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ABSTRACT 


The transport of pollutants in a cross-section of a small 
urban valley was numerically modelled using a particle-in-cel] 
technique which solved the two-dimensional diffusion equation. 
The model was based on data gathered in the North Saskatchewan 
River valley in Edmonton, Alberta during clear-sky, light-wind, 
inversion conditions. Particles, representing specified 
amounts of carbon monoxide in a roadway plume, were created 
within a V-shaped valley. These particles were advected by a 
model wind field featuring double-vortex recirculation of 
pollutants and dispersed by a simulated turbulent diffusion 
process. 

The model correctly predicted the trend to decreasing 
concentrations with increasing slope wind speed and decreasing 
source strength. Concentrations near the surface 100 m 
downslope from the source reached over 6 ppm while concentrations 
in the return flow 50 m upslope from the source reached nearly 
5 ppm. Decreasing particle density increased concentration 
variability. Concentrations were sensitive to tnitial source 
diffusion but insensitive to the magnitude of eddy diffusivities 


away from the source. 
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4.8b: Particle distribution at t=7218s produced by 


parameters of Table 4.1, except that 
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4.8c: Vertical profile of concentration 100 m 
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CHAPTER | 


INTRODUCTION 


Le legetntroduction 


Many cities in the world are located near major river 
valleys. These cities were built at a time when waterways 
were often used for travel. As the cities grew the river 
valleys formed natural cores for settlements. Today, different 
parts of the valley are often reserved as sites for recreation, 
as areas of industry, or as corridors for transportation. 
Therefore, urban river valleys have the potential for bringing 
large numbers of people together. Because of their 
micrometeorology, urban valleys also have the potential for 
trapping large concentrations of pollutants. These two factors 
could produce a localized health problem. For this reason the 
study of the micrometeorology of urban river valleys is important. 

There are few published studies concerned with small valley 
micrometeorology although some observational work does exist 


(see, for example, Klassen, 1962). More recently an experimental 
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program was undertaken at the University of Alberta in Edmonton. 
Some of the results derived from this project can be found in 
Paterson (1978), Hwang (1978) and Hage (1979). Although 
Substantial progress has been made in understanding the 
micrometeorology of the small valley, few of the new hypotheses 
have been tested. It is evident that computer simulations must 
bridge the gap between the present knowledge of valley circulation 
and future progress. 

This thesis is part of a three-way plan to model the 
unique micrometeorology of a small urban river valley. Part 
One is an attempt to model the evolution of a temperature field 
throughout the valley. This has been undertaken by di Cenzo 
(1979). In part two, Stovel (1979) attempts to show how a 
valley wind system evolves from the time-dependent thermal fields 
of part one. This thesis is part three of the plan. The author 
attempts to model the dispersal of pollutants in a small urban 


valley utilizing the valley wind system of part two. 


1.2 The Valley and its Micrometeorology 


Observations in the North Saskatchewan River valley (a 
small urban valley) in Edmonton provided data on which the 
present model is based. These data were not used in verification 


because of the model's simplicity, but rather as clues to model 


development. As a first step in modelling pollutant dispersion 
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within the valley a summary of the valley structure and known 
micrometeorology is presented. The North Saskatchewan River 
valley is typically about 1 km wide and 50 m een It meanders 
through the center of the city in a generally north-eastward 
direction. Typically, the meandering causes one wall of the 
valley to be rather sheer from the river's edge to the ridgeline. 
The other side of the valley contains a flood plain that serves 
variously as an area of recreation, industry, residence or 
transportation. It is likely that this valley asymmetry 
contributes to its unique micrometeorology. Figure 1.1 is a 
diagram of a representative cross-section of the valley, along 
with the model approximation of it. 

Evidence that the microclimate of the valley is different 
than that of the city is found in Hage (1972). Minimum 
temperatures in the valley were observed to be comparable to 
those measured at rural stations under clear sky conditions. 

The effect of the valley was to cancel the heat island effect 
of the city in this instance. 

Some understanding of the nature of circulation induced 
in small river valleys can be gained by investigating mountain 
and large valley winds. Observations showed that well-developed 
local circulations with marked diurnal variations are formed 
in valleys leading into mountain ranges. During the day the 
winds were upslope and upvalley; at night the winds were downslope 
and downvalley. The slope winds were always initiated first 


and the valley winds followed. Defant (1951) described an 
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idealized version of the diurnal change in valley winds. 
Observations in the North Saskatchewan River valley by 

Klassen (1962) and Hage (1979) generally confirmed those made 

in mountain valleys. During periods of strong prevailing winds, 

or lighter winds and overcast skies, no valley wind regime 

developed. The vertical temperature structure was essentially 

isothermal; the air in the valley was well-mixed with that of 

the plains above. When the sky was clear and the prevailing 

wind was light, the slope (drainage) wind began at about sundown. 

Initially the speed of the slope wind was small and the depth 

of flow was shallow. Both increased to their maximum values 

within two hours and persisted at or near those values until 

sunrise. The fully developed slope wind was variable. It rarely 


and its average value was about 0.2 wa 


exceeded 1.0 ms 
The depth of flow varied from about 5 m to 15 m or more. 
Observations of the temperature profiles in the valley 
revealed distinct differences between the two valley slopes 
under calm, clear sky conditions. Both slopes developed inversions, 
but that slope first receiving shading developed an earlier and 
a more intense inversion. Inversions over both slopes were 
more intense than that which developed over the urban plains 
above the valley. 
The vertical shear of the slope wind has also been 
investigated (Paterson, 1978; Hage, 1979). In well-mixed 


conditions the wind speed increased with height above the slope. 


As the stability increased, the wind speed tended to become 
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constant with height. Under the most favourable conditions, 

the speed of the slope wind decreased with height and the 

height of the maximum wind speed existed below 1.0 m. It was 
thought that this regime extended to the height of the ridgeline 
where once again the wind speed increased with height. 

Klassen (1962) found evidence of a double-vortex 
circulation system within the valley. Fog which formed after 
rain or hail drained down the slope into the river valley. 

Above the river some of the fog rose vertically and flowed 

back toward the rim of the valley. The remainder of the fog 
subsided laterally and then flowed down the valley. Therefore, 
the trajectory of a particle caught in the circulation system 
appeared to be a helix with its: axis oriented downvalley. 
Paterson and Hage (1979) estimated that the time required for 

a particle to complete one circuit of this valley vortex by 
advection alone is 40 minutes to 2.5 hours. 

Measurements of carbon monoxide (C0) concentration showed 
evidence of a distinct valley circulation (Hage, 1979). On 
most days a typical CO trace had a peak in concentration near 
1700 MST. This corresponds to the late afternoon rush hour 
traffic maximum. On some days a second peak in concentration 
existed near 2100 MST, well after the traffic maximum. At the 
downtown monitoring station the second peak was usually secondary 
in importance compared to the traffic peak; whereas, in the valley 
the second peak was dominant. The evening maximum tentatively 


can be attributed to the development of a surface-based inversion 
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formed by radiative cooling. If the transition from unstable 
to stable stratification occurs rapidly enough, the decrease 

in vertical spread of the CO plume may overcompensate for the 
decrease in source strength after the late afternoon traffic 
peak. This can result in an evening concentration maximum. 
Because the inversion is more intense inside the valley the 
evening peak in concentration will be relatively larger there. 
Most of the time the magnitude of the peaks is greater at the 
downtown station because this station is in the midst of the 
large downtown area source. Occasionally the peak concentrations 
were larger at stations inside the valley than at the downtown 
station even though source strengths are much larger downtown 
(Hage, 1979). A tentative cause of this exceptional occurrence 
is the recirculation of pollutants made possible by the valley 
winds, in association with the intensification and deepening 

of the inversion. The magnitude of the evening concentration 
peak is indirect evidence for the existence of a double-helix 
circulation regime for the valley winds. 

The object of this thesis is to investigate the effects 
of the proposed circulation system on the CO distribution 
within a small urban valley. The magnitude of the effect wil] 
be of prime importance in elucidating the proposed circulation 


regime in areas of the valley inaccessible to measurement. 
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CHAPTER I1 


THE MODELS SS THEORY, 


Dene EOcuCcLIon 


Traditionally two approaches have been used to solve 
fluid dynamics problems. In general these methods subdivide 
the region of interest into smaller cells and transform the 
differential equations involved into suitable finite-difference 
form (Gifford, 1975; Lamb et al., 1974). 

In the Lagrangian approach the grid of cells is embedded 
in the fluid and moves with it. Masses and velocities are 
usually defined at the cell corners while other time dependent 
properties of the fluid are defined at the cell center. Each 
cell corner is always associated with the same part of the fluid. 
Therefore, the cells are distorted as the fluid undergoes 
distortion. This can occur when, for example, large velocities 
exist normal to a solid boundary (Amsden, 1966). 


In the Eulerian method the mesh of cells remains fixed 
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relative to an observer while the fluid moves through the cells. 
In its strict application each cell is homogeneous withrespect 
to all fluid properties (Hanna, 1975). Cells may vary in size 
but because the cells are fixed in space the method is unable 

to resolve fine-scale structures that move with the fluid 
without introducing a fine-scale grid throughout the entire 
region. The method also produces a fictitious diffusion. 
Fictitious diffusion arises when errors in mass or concentration 
within a cell are introduced by the assumption of homogeneity 
(Lange, 1978). Since, in gradient-transport theories, diffusion 
is related to concentration, errors in diffusion rates also 
arise. This diffusion is not real but rather a by-product of 
the Eulerian method. 

In an attempt to utilize the desirable aspects of both 
methods, the present model is based on the Particle-in-Cel 
technique (PIC). The PIC method is a hybrid; it has elements 
of both the Lagrangian and Eulerian modes. The region through 
which the fluid moves is subdivided into a finite grid of Eulerian 
cells which are fixed in space and time. In general, velocities 
are defined at cell corners and concentrations are definedat 
cell centers. Masses (concentrations) and velocities are 
allowed to be time dependent. In addition the fluid itself 
is represented by Lagrangian particles or mass points which 
move through the Eulerian mesh. These particles act as markers 


for the fluid motion. Each particle is assigned a mass and the 
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total mass or concentration of a cell is determined by the 
number of particles inside the cell. For a complete description 
of the Particle-in-Cell technique, as well as its applicability 
to a wide variety of fluid dynamics problems, see the references 
Welch et al. (1965), Amsden (1966), and Lange and Sherman (1977). 
The PIC method has been applied specifically to the study of 
the dispersal of atmospheric pollutants by Lange (1973; 1978). 
A particular PIC model developed by Lange (ADPIC) has been 
assessed as a regional model in Alberta (Padro, 1979; Reid et al., 
197 2) 

Following the example of Lange, the present model solves 
the two-dimensional nonlinear transport-diffusion equation in 
its flux conservative form, relying on a given advection field 
varying in time and space, complex terrain, and time- and space- 
varying diffusion coefficients with only minor modifications 


to the computer code. 


2.2 General Description 


The nonlinear transport-diffusion equation can be 


written in the following form 
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where xy is a scalar concentration, K is the eddy diffusivity 
and U, is the given advection velocity field. Equation (2.1) 
can be simplified if the flow can be assumed incompressible. 
In the atmosphere, for small values of the Mach number, the 
only compressibility effect of importance is that related to 
the change of density with height. The domain of the model was 
the lower 50 m of the boundary layer. Therefore, the density 
can be assumed to be constant. For velocities of the order of 
l ms”! that were observed in small valleys under inversion 
situations, the Mach number will be very small. Therefore, 
the assumption of incompressibility is a good one. 


If the flow is incompressible, then 


b,ty-%- Qt 


. » (2.2) 


Combining the divergence terms the transport-diffusion equation 


(22 Wmcanmbe written. as 


— Heels a © (223) 


where Ui, is the pseudo-transport velocity and is defined as 


follows 


USM oo Uh, (2.4) 


where Up is a diffusion velocity. 
Each timestep of the model was divided into an Eulerian 
and a Lagrangian part. The Eulerian part consisted of summing 


the particles in each cell, determining the particle 
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concentration and calculating the diffusion velocity. The 
diffusion velocity was then added to the advection velocity to 
yield the pseudo-velocity. In the Lagrangian step each marker 
particle was transported along a pseudo-velocity streamline. 
Particles transported outside the grid system were counted as 
destroyed. Particles transported into the ground of thevalley 
were reflected back into the fluid proper. Based on the new 
particle positions, new Eulerian concentrations were calculated 
and the cycle was repeated. 

Cell particle concentrations were defined at cell centers 
and the velocities Uo» U, and Up were defined at cell corners. 
These velocities were then interpolated to each particle 
position. The cells in the present application were rectangles 
of uniform size. The locations of the particles representing 
fluid motion were defined by their coordinates within the fixed 
grid. The specific mass of the marker particles was assumed 
to be equal to the specific mass of air. Therefore, deposition 
was non-existent; the marker particles delineated a plume of 


pollutants. 
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2.3 Boundary Conditions 


Because of the hybrid nature of the PIC method, the 
boundary conditions were broken into two sets, one imposed on 
the Eulerian velocity field and the other on the Lagrangian 
particles. Both sets must be consistent with each other. For 
example, if there is outflow at the boundaries, then the 
particle flux must accurately represent the mass flux of the 
pseudo-velocity field. The concentration field must be smooth 
enough when it reaches the boundary so that boundary velocities 
can be specified by assuming a constant flux of particles 
through the boundary. 

The particle boundary conditions were simple. When a 
particle passed through the upper or side boundary, it was 
counted as annihilated. When a particle passed through the 
lower (valley surface) boundary it was reflected. 

Two basic boundary conditions were imposed on the 
pseudo-velocity field. They were the mass flux x u, = 0, 
corresponding to reflection of the particles at a boundary, 
and x U, = constant, corresponding to inflow and outflow of 
particles through the boundary. To be consistent with the 
particle boundary conditions, the mass flux was zero at the 


valley surface. The flux of mass was always outward at other 


boundaries. 
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2.4 The Advection Velocity 


It was mentioned earlier that observations made in the 
North Saskatchewan River valley in Edmonton and in mountain 
valleys in the United States and Europe indicated the likely 
existence of a double-vortex circulation system in the valley 
during inversion situations and that this circulation was 
largely separated from the large-scale flow. The model attempted 
to duplicate these observations by assuming that the advection 
term in the pseudo-velocity forms a double-vortex circulation. 

A typical nighttime inversion situation was modelled. 
The inversion was fully developed, that is, it filled the entire 
valley. Therefore, the double-vortex circulation also filled 
the entire valley and was assumed to be independent of the 
large-scale flow. The flow was downslope along the valley sides, 
upward (by continuity) near the valley center and then outward 
toward the valley ridge. Since the model was two-dimensional, 
the downvalley component of the observed winds was ignored. 

Several methods exist for generating valley winds. See, 
for example, Tang and Peng (1974), Tang (1976), and Stovel (1979). 
Most solve a set of differential equations including the momemtum 
equations and the heat equation in some suitable fashion. 
Considerable complexity is introduced if account is taken of 
the nonlinear coupling of the horizontal and vertical equations 


of motion and the heat equation. This complexity was thought 
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to be beyond the scope of the present model. 

Instead, a simple two-dimensional mass-conserving flow 
was sought. This can be satisfied by writing the velocity 
components in terms of the stream function wy. Then the 
magnitude of the vorticity vector ® which is everywhere normal 
to the flow is given by 


nei = M ka) 
rs xe oz] ° (255) 


For two-dimensional flow the vorticity associated with a fluid 
element is constant; in steady flow the paths of the elements 

are streamlines. Therefore, w = w(W). Provided the distribution 
of vorticity is known, solutions to (2.5) exist. This 
distribution is arbitrary for inviscid fluids (Batchelor, 1974) 
SO One may choose specific forms of wo () for which solutions 

to (2.5) are known. One convenient choice is  « wy, which 


yields the linear Helmholtz equation 


D 2 
a Se (2.6) 


Equation 2.6 is similar to that for transverse vibrations 
of an elastic membrane. Solutions are known for a number of 
shapes of membrane on which wt is constant. Of particular 
interest is the rectangular membrane which covers the domain 
G(0 < x < a, 0 < z <b). Proceeding in the normal way, assume 
a variables-separable solution, w = f(x) g(z). For the boundary 
condition » = 0, with a = 1, the eigenfunctions can be shown 


to be products of sin nmx sin maz. The complete solution is, 
a b 
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where c is an arbitrary constant and a and b are the horizontal 
and vertical dimensions of the membrane. The eigenfunction 
products form a complete orthogonal system of functions in G. 

The solution appropriate for valley winds was foundin the 
following way. Assume the valley shape to be formed by one 
half of the rectangular membrane - a triangular valley. If an 
additional boundary condition is imposed, namely w = 0 on the 
diagonal, the new solution. to (2.6) is: included: in (2.7)... tn 
other words the stream functions for a triangular valley will 
be a specific set of the complete solution. Courant and 


Hilbert (1953) show the appropriate solution to be 


yo—e tee eee, (2.8) 
NA Dal 
hates, 
Ur ec (unmrxs Ii ZT Zeer Since Sin w2)h (2.9) 
Fe) b a b 


In this case b is the depth of the valley and a is one-half the 
width. The two-dimensional non-divergent velocity field in 


cartesian components is then given by 


_~ dp 
A dz 


34 (2.10) 
A ox 


(c.6) 
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lt was a simple matter to reflect the values of the stream 
functions resulting from (2.9) about the valley axis so that 
the desired valley winds were produced. 

A comment about the realism of the velocity field as 
generated by (2.9) and (2.10) is in order. The field was 
experimentally reproduced on a finely spaced grid and the 
maximum values of Un and W, were found. The ratio of the 
maximums (u,/wa) was found to be approximately 6.5. For synoptic- 
scale flow, this ratio is the order of 1000 (Holten, 1972). 

As it is in synoptic-scale flow, the ratio of horizontal to 
vertical wind speed in a valley is governed by the dimensions 

of the boundary. For example, the ratio in synoptic-scale flows 
depends on the height of the tropopause and the horizontal 

size of synoptic disturbances. A similar relationship was assumed 
to exist for the flow in a river valley. The North Saskatchewan 
River valley is approximately 800 m wide and 50 m deep. The 
expected ratio in the valley is, therefore, nearer to 16. The 
ratio in the generated velocity field can be shown analytically 
to be a/b, twice the ratio of the valley dimensions. Thus (2.9) 
yields values of vertical velocity that are larger by a factor 

of two over those expected. The finite-difference forms of 
(2.10) further enhance the vertical velocity relative to the 
horizontal. The vertical dimension of a pollutant plume 


trapped in such a velocity field should, therefore, be increased. 
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A more serious departure of the model winds from inferred 
small valley winds was the relatively large maximum in the 
horizontal velocity that occurred near the upper boundary. Model 
horizontal velocities here were comparable to or larger than 
those on the valley slope; whereas, actual velocities here are 
believed to be smaller. Thus, a plume traversing a circuit of 
the model double-vortex circulation system should undergo 
relatively less vertical diffusion. Streamlines of arbitrary 


magnitude are shown in Figure 2.1. 


2.5 The Diffusion Velocity 


All ramifications of gradient-transport theory depend 
ultimately on the notion that the flux of a quantity is 
proportional to its gradient. That there seems to be no 
precise physical basis for this assumption is emphasized in 
Slade (1968). The theory does, however, provide useful results 
and will be used to define the diffusion velocity as that 
velocity (of a pollutant particle) resulting from a concentration 


gradient. The diffusion velocity is given by 


U -- & (Ome) 


where K is the eddy diffusivity and x the scalar pollutant 


concentration. This definition follows directly from (2.4). 
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The concentration x is simply the total number of 
Lagrangian marker particles in a unit volume. (The length of 
the downstream component is arbitrary.) The finite-difference 
form of (2.11) will be discussed in Section 3.5. The remainder 
of this section will be used to discuss particular diffusion 
parameters. 

In principle the model can accomodate the full two- 
dimensional cartesian eddy diffusivity tensor ype In practice, 
however, the turbulence was assumed to be isotropic so that only 
the diagonal elements, ve and Kee survived. It then remained 
to determine forms of the vector eddy diffusivities that were 
both realistic and tractable. To do this, it is necessary to 
review some of the wind-profile observations made in the North 
Saskatchewan River valley. 

The time of interest was nighttime, after approximately 
2200 MST. The inversion was entrenched throughout the entire 
depth of the valley, as were the valley winds. The maximum 
downslope wind was found to exist below a height of 0.8 m 
(Hage, 1979). Therefore, the downslope wind decreased with 
height and, except perhaps for the lowest few centimeters, a 
constant momemtum flux layer did not exist in the valley. Indeed, 
because of the intense inversion, the valley air may be more 
accurately characterized as a zero vertical turbulent flux layer. 
Vertical velocities at this time were found to be below the 


resolution of sensitive propeller anemometers; they were 
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essentially zero. Horizontal velocities were thought to be 

much larger because of the large horizontal temperature gradient. 
These observations served to put severe limitations on the 
functional forms of the eddy diffusivities. 

In an attempt to put bounds on the eddy diffusivities, a 
range of values were used (see Section 3.5). The vertical 
coefficient ranged from near molecular to near that predicted 
for a stably-stratified atmosphere assuming constant flux. 

The horizontal coefficient ranged from several orders of 
magnitude larger than molecular to a value consistent with 
stable stratification. For all cases except those assuming a 
constant flux, directional variations of the eddy diffusivities 
are not known. Therefore, constant values were used. 

Functional forms of the eddy diffusivities have been 
postulated assuming the flux in the lower 50 m to 100 mto be 
constant. Even though this assumption is inferred to be 
incorrect for small valleys, the model retains the capability 
of computing eddy diffusivities from it. In the determination 
of a form for the vertical diffusivity, turbulent diffusion 
induced both by drag on the earth's surface by the atmosphere, 
and by thermal convection within the atmosphere must be taken 
into account. Assuming the inertial subrange hypothesis to be 
valid for the region of the boundary layer of this application, 
similarity theory can be used to determine the vertical diffusivity. 


A tacit assumption is that the flat plain forms of the eddy 
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diffusivities are approximately valid for a shallow to moderately 


sloped valley. With this approach a can be shown to be 


= Uk é 
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K 


na (271 2) 
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where U, is the friction velocity and k is the von Karman 


coefficient. The dimensionless wind shear o is given by 


(Zea) 


where U is the average horizontal wind speed. Semi-empirical 
relations have been used to determine values for bat Businger 
(1973) gives the following relations based on the AFCRL Kansas 
experiments: 
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where t is a measure of atmospheric static stability and is 


defined by the relation 
c= 2/L (215) 


where z is the vertical height and L the Monin-Obukov length. 

In order to be internally consistent one should calculate 
both U and tc for use in determining Ke: However, the calculation 
of t was thought to be beyond the scope of the model. In 
practice a value of static stability was specified asmodel input 


and U,was calculated (and used as model input) from gathered 
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data. 

Functional forms of the eddy diffusivity can be found by 
making assumptions about the shape of the wind profile and the 
vertical variation of the vertical eddy diffusivity 


(Pasquill, 1974). For a power-law profile, that is, 


2A) 


where K and U, are values at a fixed reference height Z)> 


Davies (1950) gives the relationship 


K Zier Oxi.) 


lt has been observed in an atmosphere of neutral stability 
with a power-law wind profile, the value of the exponent m is 
approximately 1/7. With n = I-m (the conjugate power-law 
condition), (2.17) results in a ground-level axial concentration 
varying as aloes whereas, the observed dependence is ee 


for a continuous point source at ground level (Pasquill, 1974). 


: . a 
To obtain this observed dependence, Kon must vary as z , where 


a = 1 + m(1-3m)/(l+m) . (23113) 


For m = 1/10, indicative of a slightly stable atmosphere, a is 


approximately 1.07. 
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As a first approximation to the eddy diffusivities ina 
stably-stratified atmosphere, it might be assumed that eee and 
Nee are of the same order of magnitude. If, as was supposed 
ear liver seUneand on are constants, then Sop « EEL With this 
vertical dependency, oe and Se can be of the same order of 


magnitude for all z within the range of the model. 


Equations (2.12) and (2.17) have an advantage inthat 
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neither is directly dependent upon the plume parameters or and O,- 


This is valuable when marker particles from different sources 
(and in general different plume dimensions) intermingle. To 
treat them as separate plumes would be unrealistic in this 


situation. 


2.6 Sub-grid Scale Diffusion 


It can be shown (Lange, 1973) that when the particle 
distribution cannot be resolved by the grid mesh, errors in 
diffusion rates occur. This happens when continuous sources 
much smaller than the grid spacing are modelled. The problem 
is compounded in the presence of advection because the errors 
of inadequate resolution are carried downstream. 

This problem was treated in the following way. Theinitial 
shape of the plume was assumed to be Gaussian. The particles 
were dispersed, the plume retaining its initial shape, until 


the plume could be resolved by the grid mesh. At this point 
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the locally-Gaussian assumption was discarded and diffusion was 
determined by the grid point concentration gradient method. 

The Gaussian diffusion velocity can be found in the 
following way. Assume the concentration field in one dimension 


to be given by 


-X 2 
own | 
xX =e exp (2.19) 
2710 aos, 


Where Q is the source strength, om the plume standard deviation 
and Xp the distance from the particle to the plume center of 


mass. From the definition of the diffusion velocity 
"Ki d2n 
Ws ce Oe eS, (2411) 
D Bee GbR 


Substituting from (2.19), 


p= a 03 (2.20) 


rife pane (2221) 


Equations (2.20) and (2.21) can be combined and integrated 
over the length of the timestep to yield a diffusion distance 
AXy- However, assumptions must be made about the time- and 


space-dependence of Net and OL. For a constant diffusivity, 


K = constant, 
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where omy is the initial standard deviation in the x-direction. 


For this case the diffusion distance can be shown to be 


re At Lise 
Bin iodine all bc yeenera creme reg rp (2.23) 
Ox XX 


For three-dimensional, scale-dependent diffusion, Walton (1973) 


gives the following relations 


Peete oa (2.24) 


Xx x 
o 2 = (o 2/3 = 2/3ce!/ 34)? (25253) 


where ce is the rate of eddy energy dissipation and c is a 
constant of the order of 1. Because the vertical turbulence 
is to a large extent suppressed, a case can be made for using 
Lin's (1972) two-dimensional forms of (2.24) and (2.25). 
However, they were not used here. With c = 1, the diffusion 


distance can be shown to be 


AX, = Xp ih oe ONE Bie 
Oo 1 OF (2.26) 
1.5 ox}1/3 — 
‘ € 


Similar relations hold for the vertical direction. The total 


distance travelled by the particle is then 


AX = AX a unAt. ee) 


It should be noted that the sub-grid scale diffusion scheme was 


used only to produce a realistic source configuration on a 
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scale that can be resolved by the grid mesh. The grid was 
assumed to resolve the plume when the standard deviation of 


the plume becomes greater than one cell length (i.e. o, > Ax). 


Xx 


2.7. the source 


The source under consideration in this model was a plume 
resulting from automobile exhausts. Although hydrocarbons, 
particulates, carbon dioxide, water vapor, and oxides of sulphur 
and nitrogen are also exhaust components, only carbon monoxide 
was considered. In addition it was assumed that CO is a 
chemically inert species, although this is not precisely true 
(Sandhu, 1975). 

In general the diffusion of pollutants from motor vehicles 
depends on the following factors (Fanaki and Kovalick, 1974): 

1) number, type and age of the vehicle used, 

2) geometry and configuration of roadways, 

2)memissionumate (aetUnction ofevenicle speed), 

4) vehicle aerodynamics and vehicle spacing on 

the road, and 

5) atmospheric variables. 
Included in 5) above are factors such as wind speed and direction 
and atmospheric stability. Factors 1) and 3) can be approximated 


by using an appropriate value of the emission factor, as in 
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Taylor (1973). By restricting the model to two dimensions, the 
complexity of factor 2) is reduced. See Section 3.6 for details 
of the roadway configuration. 

The last two factors contribute to the shape of the 
exhaust plume in two ways. Excess buoyancy is created by a high 
exhaust temperature (a typical value is 227 C). The height of 
rise, and so the vertical dispersion of the plume, are dependent 
on stability. One means of dealing with this is by determining 
an effective source height, where the effective height is 
defined to be that height where buoyant rise can be neglected. 
In addition, mechanical mixing is important in the wake of the 
moving vehicle. Because horizontal dispersion is augmented 
by mechanical mixing, Fanaki and Kovalick (1974) assume a 
virtual point source to exist upwind of the roadway under the 
@frect Olan) across-road wind. Danard (19/2) Uses a value of 
20 ae for Se in the lowest 3 m as a means of dealing with 
increased vertical dispersion. 

In practice the present model treated factors 4) and 5) 
in a common fashion. Buoyancy was ignored and the source was 
assumed to be at ground level. The effects of mechanical 
mixing were dealt with by assuming constant values for the 


initial horizontal and vertical standard deviation of the plume. 


See Section 3.7 for the particular values used. 
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CHAPTER III 


THE MODEL - PARAMETERIZATION 


3.1 Introduction 


Chapter I! contains the theory that forms the basis of 
the model. In this chapter, model parameters are introduced, 
finite-difference forms for the equations are inspected and the 


calculations used in the computer code are reviewed. 


3.2 The Valley 


Although the model can, with minor modifications to the 
computer code, accomodate a valley of any well-behaved shape, 
the valley shape chosen was a simple V. This simplified the 
production of an advection velocity featuring return flow and 
was thought to be a reasonable starting point in determining 


the suitability of the Particle-in-Cell approach for modelling 
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30 
valley winds on this scale. The model used a valley depth of 
50 m and a width of 800 m although they were input variables 
and could be changed to suit the application. These dimensions 
are approximately those of the North Saskatchewan River valley 
in Edmonton. Experiments at this site provided both input data 
and verification data. Summaries of the experimental results 
can be found in Paterson (1978) and Hage (1979). 

Because the valley wind system was assumed to fill the 
valley, the model assumed no large scale flow; in this 
application the model was independent of the overlying flow. 
Thus, the valley may be oriented in any direction. The domain 
of the model extended from the valley bottom to the ridge in the 
vertical and from ridge to ridge in the horizontal. Again, the 
model has the capability of an extended range in the vertical 
by altering the appropriate input parameter. This will 
necessitate the specification of an overlying flow. 

The choice of grid size depends on the scale of the flow 
to be modelled and on the scale of the turbulence to be 
resolved. The grid size (along with the length of the timestep) 
was also chosen so as to minimize efficiently the truncation 
errors inherent in the finite-difference schemes that were used. 
The cells were chosen to be of uniform size throughout the grid. 
They were rectangles and varied in height from 2 to 5 m and in 
length from 20 to 40 m, although the size was constant during 


each run. Again, grid size was an input variable. 
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343 “Cell Areas 


Subroutine SRAREA of the program (Appendix B) computes the 
sizes of all cells in the grid. This allows the model to use a 
valley of arbitrary shape. Associated with an arbitrarily shaped 
valley are cells of arbitrary shape. The code assumes, however, 
that the valley surface is linear between grid points. Since 
the valley was chosen to be V-shaped, minor modifications must 
be made to the computer code to allow for an arbitrary valley 
shape. In particular, the slope which in this application 
is constant must be allowed to vary. 

The cell mesh was set up throughout the entire domain of 
the model. The cell number was chosen to be a two-dimensional 
array with origin at the bottom left cell coordinate. In this 
way cell number was linked directly to cell position. The cell 
area was defined to be that area within the valley ''atmosphere", 
excluding any part below the valley groundline. The four corners 
of each cell were checked to see whether they were above or 
below the valley surface. The parameter N was introduced to 
facilitate this. Beginning at the lower left corner, N was 
increased appropriately if the corner in question was above the 
valley surface. See Figure 3.1. For example, if the upper 
right corner was in the free atmosphere, N = N+100. Fora 
cell totally above ground, N = 1111. Based on the value of N, 
an appropriate scheme was used to calculate the cell area. 


The area of a ''complete'' cell was simply the product of horizontal 
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and vertical grid spacing. The area of a cell totally below 
groundline was set to -1. The area of ''incomplete'' cells was 
found by finding the intersection of the cell walls and the 
valley groundline. Since the valley was linear between grid- 
points, the area was found by the appropriate summing of the 


areas of smaller rectangles and triangles. 


3.4 The Advection Velocity 


According to the experimental results in Paterson (1978), 
an average value for the magnitude of the downslope wind in 
the time period from 2200 to 2400 MST was approximately 0.2 on 
although the real wind is quite variable. By varying the 
constant in (2.9), it was possible to produce advection velocities 
of this magnitude. The difference between downslope and 
horizontal was neglected because of the shallow slope in this 
application. Setting the constant C = 4.0, the maximum value 
of Ju, | approached 0.5 ws wa | approached 0.07 ms | This 
ensured a reasonable value for Uns Note that although the 
wind field in the present application was steady state, 
considerable complexity could be introduced by choosing the 
"constant'! C to be time-dependent. In fact, Defant's (1951) 
diurnal variation of an idealized valley wind could be approximated 
in two dimensions. If the height of the inversion and the eddy 


diffusivities are also allowed to be time varying, then the 
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valley wind system can be made quite realistic. 

The vertical shear of the horizontal advection velocity 
was also investigated for comparison with that observed in the 
North Saskatchewan River valley. Figure 3.2 shows the results 
for a position near the middle of the valley slope. The 
advection velocity assumed by the model had a vertical shear 
that was similar to the observed shear for the time period in 
question. A similar shear existed at the upper boundary even 
though this may not be realistic. The model wind had its 
maximum velocity at the valley surface; whereas, the real 


wind has its maximum below 0.8 m but above a height of several 


centimeters above the surface (the viscous sub-layer). In this 


respect, too, the model wind was not entirely realistic. 


As mentioned in Section 2.2, the advection and diffusion 
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velocities were defined at cell corners. The advection velocities 


were derived from the stream function field using central 


difference approximations 
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where Ax and Az are the cell dimensions in the horizontal and 
vertical. The truncation error of uy (w.), that is, the 
magnitude of the remaining terms in a power series expansion 


of w(x,z), is of the order of Ax2 (Az7) (Haltiner, 1971). In 
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Figure 3.1: Parameter N used in determining 
cell areas. 
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Figure 3.2: Experimental and model wind profiles. Solid 
lines are typical experimental profiles (from 
Hage, 1979). Dashed line is model wind profile. 
Local time. 
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general, the higher the order of the truncation error, the 
more accurate is the finite-difference approximation. 
The program is capable of producing any desired wind 
speed by the choice of a value for the constant C. The user 
may also specify the advection velocities to be everywhere 
zero; that is, dispersion will take place by turbulent diffusion 


only. 


3.5 The Diffusion Velocity 


The diffusion velocity was defined in Section 2.4as 


ee (Onl) 
A range of values for K are considered. Determine first an 
upper bound for the eddy diffusivities. On a flat plain ina 
stably stratified atmosphere (2.12) and (2.17) are valid. 
However, the presence of the valley adds complications. Two 
factors tend to limit diffusion. Firstly, the dimensions of 
the plume are controlled ultimately by the size of the valley. 
Secondly, the rates of horizontal and vertical diffusion are 
decreased by increased stability (Tennekes and Lumley, 1972). 
Since the inversion is more intense in the valley, diffusion 


there should be retarded relative to the flat plains. The 


relative magnitudes of these two factors in a small valley are 


‘uncertain. However, for simplicity it was assumed that they 
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cancel approximately any increase in dispersion which may be 
due to complex terrain in the valley. Therefore, as an upper 
bound consider the case of a slightly stable atmosphere in a 
constant-flux layer. Then (2.12) and (2.17) are approximately 
valid. As a first approximation, it is desired that Ke and 


KS be of the same order of magnitude. Therefore, assume 


LL pe (3.3) 
K = kK ° = 
XX ZZ We0 ) 

7h mM 


The value of the von Karman coefficient k was set to 0.35. This 
procedure (3.3) is valid only for a time period near 1800 MST 
(see Figure 3.2), and was used only to find an approximate upper 
bound. A value of 0.1 e was suggested for U, although it is 
an input variable. Because the non-dimensional shear tn isa 
function of stability ¢, the computer program also retains f 

as an input variable. A value of ct = 1.05 was used in this 
application. Use of these particular numbers resulted in values 


Oiakempancak@erOretherorder on Onl feel However, in the 
XxX Zz 


time period under consideration the valley cannot be realistically 


thought of as having a constant-flux layer, as these equations 
imply. Without this simplifying assumption, formulations for 


Sele and Kee are probably more complex and are at present 


unknown. It has been shown, however, that an upper bound for 


2 iz A further 


the diffusivities should be of the order of 0.1 m°s- 
clue to their magnitude is that a large horizontal temperature 


gradient exists near the slope. With this in mind, an upper 
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bound for the diffusivities is suggested to be ors =eOel mis) 


and Se mS ORX lo? eet These were assumed to be constant 
throughout the valley. It should be stressed that, although 
the code was capable of calculating yd and Soe directly from 
(ala )ymande( 2.7 jemchisatechni quem!s valideonly vonsastiateplain. 
Because forms of the diffusivities are unknown for small urban 
valleys, the present model experimented with a range of constant 
values. 

Because the observed vertical velocity was very small, a 
reasonable lower bound for See may be the kinematic viscosity 
vy. At a temperature of 10 C and a pressure of one atmosphere 
v is approximately 1.8 x Os | 

A lower bound for idee was difficult to establish because 
guidelines do not exist as they did in the case of the upper 
bound. However, a much higher value than that for oF was 
expected because of the large horizontal temperature gradient 
near the slope. Using this fact as a guide, the model assumed 
a lower bound of Kes al. Oinex eee 

As stated earlier diffusion velocities were defined at 
cell corners and concentrations at cell centers. Therefore, 


x was defined as the average concentration of the four cells 


adjacent to each corner. In the notation of Figure 3.3 


X= W(x) + xy + x3 + Xy) (3.4) 
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Consider the one-dimensional gradient MEX" In finite- 
difference form the gradient is written as 


oy ag CG = Dep se X3) (3.5) 
2Ax 


is 


Therefore, the one-dimensional diffusion velocity Un. 


written as 


i 2K x (x5 te X3 7% xX] = Xy,) 


Dx Ax (x, POG) + X3 + x)) (3.6) 


Similar expressions are valid for VEX and U Equation 3.5 


Dz° 
is accurate to second order in Ax. 

The maximum diffusion velocity occurs when a single 
particle is surrounded by empty cells. The maximum velocity 


allowed by the computer code was 


Uae see? Can 


Ax 

A feel for the magnitude of the diffusion velocity can be 
found with the following example. For a horizontal eddy 
diffusivity of magnitude 0.01 ia and a horizontal grid 
spacing of 20 m, the maximum diffusion velocity was approximately 
| mm a 

The diffusion velocity for cell corners below the valley 
groundline was set to zero. This acted to decelerate artificially 


particles diffusing toward the ground. 
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The advection velocities and diffusion velocities (when 
the particle was not considered part of a sub-grid scale plume) 
were interpolated to the particle positions from neighbouring 
grid points using a two-dimensional linear interpolation scheme. 
The forms for the horizontal u and vertical v velocities are 


given below in the notation of Figure 3.4: 


nes x, (U52Z, + U,Z,) + x, (U,2, + UZ) ey 
AxAz : 


Z,(WiX, + WX) + z, (W x, + W x1) 


aa ay | he) 


3.6 The Source 


It was stated in Section 2.7 that mechanical mixing due 
to vehicle motion is treated by assuming an initial source 
standard deviation. Zimmerman and Thompson (1975) assumed the 
standard deviation of the plume in the horizontal to be 
approximately equal to one-half of a car length and that the 
vertical standard deviation was somewhat less. The same logic 
was followed in this application. The initial horizontal 
standard deviation of the automobile plume was set to 3 m and 
the vertical to 2 m. Again, both parameters were input variables. 


The code used a Gaussian random number generator to produce 
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Figure 3.3: 
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Figure 3.4: 


Scheme for calculation of 
diffusion velocity Up: 
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Velocity interpolation scheme. 
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particles of this configuration. These plume dimensions were 
valid only for the initial particle placement at each timestep. 
The variance of the plume at later times was calculated as the 
average of the sum of squared particle distances from the 
center of mass of the plume. 

The number of particles produced at each timestep depended 
upon the traffic count at each source location. Ultimately 
the number depended upon the storage capacity of the computer 
and the expense of running the program. In its present form 
the code produced, On average, about | particle each second of 
model time. This was about 7000 particles for a two hour run. 
Naturally, the more particles produced at each timestep the 
better the resolution. Concentration and particle position 
fields were smoothed with increasing particle numbers. However, 
the time and cost of running the model with large numbers of 
particles can become prohibitive. 

Each source in the model was essentially a point source 
(actually an area source because the plume had a finite initial 
dimension). Source coordinates were input to the program as 
variables. All sources at present were assumed to be on the 
surface of the slope. Since the model was two-dimensional, 
each point source represented an infinite line source running 
the length of the valley. 

The mass per unit time of CO produced at each point source 
is given by 
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42 
where EF is the vehicle emission factor and D is the length of 
a segment of an infinite line source over which the emission 
factor is valid, NV is the total daily traffic count along the 
roadway, F is a fraction accounting for the diurnal traffic 
cycle and G is the mass of CO produced per unit time in km. 
Paterson (1978) used a value of 23 g tay per vehicle for the 
emission factor in approximating the traffic mix in the North 
Saskatchewan River valley. The same value was used here. The 
distance D was arbitrary but necessary when converting from 
model concentrations (particles per cell) to commonly used 
concentration units (ppm or a. A value of 50 m was assumed 
by the model. 

The product NV*F yields the total number of vehicles 
passing through the distance D in one second. In this application 
the plume of only one source was modelled although the code set 
to ten the limit on the number of point sources that may be 
input. The point source modelled a line source such as River 
Road which runs parallel to and inside of the valley for a 
distance of several kilometers. The City of Edmonton (1978) 
provided traffic count data for River Road. These data were 
used to determine values for NV and F. The specific value used 
for NV was 16400 vehicles/day, a typical weekday total. A 
linear variation of traffic count with time was assumed for the 


hours 2200 to 2400 MST. With reference to Figure 3.5, the 


fraction & 1S Of. form 


_ 0.04 - 0.0075T/3600 


3600 (3.11) 
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where T is model time in seconds. Hour zero model time is 
2200 MST. The units of F are s_! 

[If one Lagrangian marker particle is equivalent to 
0.20 grams of CO, then for a distance D of 50 m and cell 
dimensions as previously specified, the number of Lagrangian 
particles produced at each timestep is given by 

N = 26.19 At (0.04 - 0.0075T/3600) (3902) 


where At is the length of the timestep. 


3.7 The Timestep 


The choice of a timestep length is dependent upon several 
factors. Minimizing the timestep will minimize the truncation 
error inherent in the finite-difference schemes and will 
provide for smooth fields of concentration and particle position. 
Increasing the timestep will decrease computer costs. A third 
factor is also involved. Often one wishes to resolve ''waves'' 
or disturbances of a certain size in the concentration field. 

In this case, the length of the timestep should depend upon the 
"wavelength'' to be resolved and on the particle velocity. 
The model used an approach of this sort. The timestep was 


a function of the grid spacing and the maximum velocity components 
(Lange, 1973) 
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This ensured that the maximum distance travelled by any particle 
was one-half of a grid spacing in any direction. Since the 
concentration field was found by summing the particles ina cell, 
the maximum resolution in the concentration field was the Nyquist 
wavelength - 2 grid lengths. The velocities ite and Wren 
the maximum velocities experienced by a particle. They are not 
velocities at grid points. The code allowed the timestep to 
be decreased as necessary from step to step but allowed an 
increase in the timestep of only 30%. This was to aid in 
stabilizing the model (Lange, 1973) although as a precaution 
only. The problem was not encountered here. It should be 
noted that the length of the first timestep must be input; the 
henotn ofeali@later steps is calculatedras in (3-13). 

Once the length of the timestep had been calculated, the 
particle was advected along a pseudo-velocity streamline. The 


scheme used to do this was a forward timestep of the form 
XP(t + At) = XP(t) + us At (3.14) 


where XP is the x-coordinate of the particle. The scheme is 
first order with respect to At. The calculation of a new 


vertical coordinate was identical. 


3.8 Particle Trajectories and Reflections 


Where a particle was transported below the groundline, 


the code reflected the particle back into the valley atmosphere. 
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45 
In order to accomplish this the intersection of the particle 
trajectory and the groundline must be calculated. Thiswas 
done in the following way. The slope of the trajectory could 
foe calculated since the coordinates of the beginning and end 
points of travel were known. The slope of the groundline was 
given. Using an arbitrary point on each of these lines, the 
z-intercepts were determined and so the equation of each line. 
The calculation of the coordinates of intersection was then 
straightforward. 

For a valley shape other than a simple V, the local slope 
must be estimated. In this case it was assumed that the 
distance travelled by the particle was small. Then the 
groundline between beginning and end points of travel could 
be assumed linear. 

Once the intersection was known, the final, above-ground 
particle position could be determined. Refer to Figure 3.6. 

The distance AB could be found from the beginning and end points 
of the particle motion. The coordinates of the trajectory- 
groundline intersection were known. Therefore, the lengths OA 

and OB could be calculated. The slopes of the particle trajectory 
and groundline were known so that the angle of incidence i could 
be found. The point C was then determined, assuming the angles 

of incidence and reflection to be equal. The distances OB and 

OC were assumed equal; the ''collision'' of the particle with 


the slope was assumed elastic. 
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Figure 3.5: Weekday mean traffic. City of Edmonton: 1977. 


Figure 3.6: Reflection of particles at groundline (solid 
line). A denotes initial particle position. 
C denotes final position. 
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The above calculations necessarily involved the 

determination of the position of a particle relative to the 
valley surface. Either a particle was below the ground]ine or 
it was not. This entailed comparison of two numbers to the 
limits of their accuracy where errors due to computer round- 
off were important. Experience showed that those calculations 
compromised the accuracy of the last significant figure. 

The code, therefore, set an arbitrary accuracy limit. In 
cases where the position of a particle relative to the 
groundline was determined incorrectly beyond this limit, the 
particle involved was simply destroyed. Incorrect determination 
of particle position within the accuracy limit resulted in 


immediate termination of the computer run. 


Oe, ee 


len imtatab *jat%Goni. ~.ba “vos seat Vigne 2ew Bey Havnj ofaisyeq. 


aig viene fj dimascn 

sty of avilelsy pis! Fiac s Te inate oat 
ya so bavore' eid aolhe faw sfoli<sa Ss agatha 
anda of Siadoun owy Tonos !séanos hehibins 21A7, 
-bavo > Isiligmos of subleto vis siotw Jahuose thot i 
srorislyupten $20n3 809 DSswore 32193 3aqKe nese 90 ts 
.890pl? 27851 Haste teal sly 3 apiOsa6 siz beeime qf “7 


a 


ol. <dimit VYoswons vieriidie née sae sD tS TANS ~2too “dT 
; carn ¢ 
SOS 9I,.OViJ Fe) SIs sg 169 8 lo acitizgog Sno acti 28862 — 


oH3,yJTinid Zins broysd *y{iS<*7o9Nt “bsnl mszeb, eeM sti ibrverg: i 


badlurei ttnl! yo usse sig Aletsiv noha fay ofai sine ms 


4 Yesuanda 343 To noi aentimtad mn a) 
+ wl 


48 


CHAPTER IV 


RESULTS 


4.1 General 


In this chapter, the sensitivity of the model is appraised 
with regard to changes in grid spacing, magnitudes of the 
vertical and horizontal eddy diffusivities, initial sizes of 
the automobile exhaust plume, magnitudes of the advection 
velocity, and particle density. 

Two methods were chosen to display predicted changes in 
CO concentration that result from variations of the above 
parameters. The first is a strikingly vivid scheme of plume 
visualization. Since the model represents pollutants as quantized 
particles, the locations of all particles within the grid can 
be displayed on a scatter diagram. Appendix A contains a series 
of these plume snapshots. Although they are very useful for 
displaying plume motion and relative density, they are of limited 


use in quantitative display. To this end a second form of 


display was used. Two points on the valley slope were chosen, 
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4g 
one 100 m downslope from the CO source and the other 50 m 
upslope from the source. At each of these positions, the CO 
concentration in parts per million of the three grid cells 
immediately above the valley surface was extracted at each 
timestep of the model. In practice, the length of the timestep 
ranged approximately from 12 s to 65 s depending on the particular 
parameters used, but was nearly constant during each simulation. 
For this reason, and to smooth the graphical output, the 
concentration in each of the chosen cells was averaged over 
approximately a three minute period. The result is a vertical 
profile of pollutant concentration as a function of model 
time at each location. 

At this point it might be asked why the entire concentration 
field is not displayed and contoured, since this field was 
calculated and utilized as part of the computer code. The answer 
is that, in general, the boundaries of the cells and the valley 
surface were not aligned, nor was the valley surface coincident 
with the cell diagonals. Since the concentration of all cells 
totally below the valley surface was defined to be zero, computer 
contouring of the concentration field resulted in an 
unrepresentative display. 

Table 4.1 contains a list of standard model input parameters. 
While some remained unchanged from run to run, others were 
allowed to vary as part of sensitivity analyses. The particular 


values appearing in Table 4.1 are those used by the computer 
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Table 4.1: Standard parameters of the computer model. 


Parameters marked 


ate 
a 


of sensitivity analyses. 


Parameter Symbol 


L 
HPRIME 
INCX* 
INCZ* 


SDEVX* 


SDEVZ* 


DKMX* 
DKMZ* 


AVWF% 


TMAX 


Description 


valley width 

valley depth 
horizontal grid size 
vertical grid size 
initial horizontal 
standard deviation of 


exhaust plume 


initial vertical standard 
deviation of exhaust plume 


horizontal eddy diffusivity 
vertical eddy diffusivity 


maximum advection velocity 


maximum length of 
simulation 


number of sources 


horizontal coordinate of 
source 


representative weight of 
each particle 


are allowed to vary as part 


Value 
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program in its standard simulation. Most sensitivity analyses 
were compared against this standard. It should be emphasized 
again that, although the computer code permits other alternatives, 
eddy diffusivities were held constant during each run. For the 
time period involved diffusivities were assumed to be independent 
of time and space. In all cases, source strength decreased 
linearly with time as in (3.11). Appendix B contains 
the complete computer program. 

The remainder of this section contains a discussion of 
two important model assumptions. The first is conservation of 
mass. Although the advection field conserved mass, the presence 
of diffusion added complications. Were the two processes 
handled separately throughout the entire code, all particles 
diffusing through the upper boundary could have been annihilated. 
This mode would have increased the cost and decreased the 
accuracy of the computations. Instead, the model handled the 
two processes, as much as possible, as one and created the 
problem of determining whether an annihilated particle was 
diffused or advected through the upper boundary. The code 
dealt with this as follows. The ratio of maximum diffusion 
velocity to maximum advection velocity was computed. For the 
parameters of Table 4.1, this ratio was typically of the order 


3 


of 10 7 to 10? or smaller. Thus, for example, every thousandth 
particle that was transported through the upper boundary was 


assumed to have been diffused there and was annihilated. Other 
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particles transported through the upper boundary were assumed 
to have done so because of the finite-difference representation 
of the advection velocity. For a two hour simulation, the 
number of particles annihilated as a result of this method 
ranged from zero to three. The total number of incidents of 
passage through the upper boundary was typically two to five 
thousand. Thus, the conservation of particle mass approximated 
the conservation of pseudo-velocity mass flux. 

The second assumption concerns fictitious diffusion. 
Lange (1973) stated that the PIC technique on which the present 
model is based eliminates fictitious diffusion which is inherent 
in any Eulerian method. However, the term "fictitious diffusion" 
was not precisely defined and this caused some confusion. In 
its basic form, fictitious diffusion is caused by computer 
round-off. Multiple calculations increase the magnitude of the 
error. This type of error is impossible to eliminate ina 
numerical model. Fictitious diffusion is also created when 
concentrations and diffusion velocities are defined in terms 
of a grid. Theoretically, by choosing the grid and the timestep 
to be very small, these errors can be reduced to near the 
magnitude of the round-off error. In this application, predicted 
concentrations were sensitive to changes in cell size (see 
Section 4.3). Other processes contributing to fictitious 


diffusion were thought to be dominated by the effects of changes 
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4.2 Gaussian and Concentration-Gradient Diffusion 


As discussed in Section 2.4, the model relied ona 
concentration-gradient technique to simulate diffusion. For 
sub-grid size particle distributions, a Gaussian technique was 
used. This section briefly compares the dispersion produced 
by these two schemes. 

For purposes of this comparison only, the source was 
positioned at the center of the valley (400 m, 25 m). The 
duration of the simulation was 20 minutes. Diagnostics indicated 
that, approximately 11 minutes after start-up, particles were 
reflected from the upper boundary. Since the source was at the 
valley center, reflection from the valley surface was expected 
at about the same time. Therefore, boundary effects will be 
minimized by examining dispersion within the first 11 minutes. 

Figures 4.1la and 4.1b show dispersion estimated respectively 
by Gaussian and concentration-gradient diffusion at approximately 
1] minutes after start-up. In both cases, vertical and 
horizontal diffusivities were set to a constant value of 


0.1 ane 


Therefore, distortion of the plume from circular was 
produced by scaling. Horizontal grid spacing was 8 m; vertical 
grid spacing was 5m. Initial size of the plume was large enough 
so that, for the case of concentration-gradient diffusion, the 
initially-Gaussian specification was unnecessary. 


Some differences are evident between Figures 4.la and 4.1b. 


The plume of Figure 4.1b appears smoother and more disperse. 
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Figure 4.1a: Gaussian dispersion of plume at t=646s. 
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Figure 4.1b: Concentration-gradient dispersion of plume at t=664s. 
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In both cases approximately 650 particles were produced (about 

|] particle each model second). However, due to the differences 
in method, the length of the timestep (and the number of 
particles produced each timestep) was larger by approximately 
20% in the case of concentration-gradient diffusion. This is 
Opposite to what might be expected by examining Figures 4.la 

and 4.1b. The method with the larger timestep might be 

expected to produce a larger particle density near the center 

of the plume. Figure 4.1lc is a plot of concentration vs. time 
of one model run for a cell with center at (428 m, 27 m). 
Plotted is the unsmoothed concentration at every fourth timestep. 
The concentration quantum was 0.025 particles/m-. Figure 4.1c 
illustrates qualitatively that the rates of diffusion produced 
by the two methods are approximately equal on the large scale. 
This suggests that when the calculation of distance diffused 

by particles transfers from Gaussian to concentration-gradient 
technique, the change will not result in a discontinuity in 

the rate of diffusion. Differences in the appearance of 

Figures 4.la and 4.1b may also be attributed to small variations 


in the rate of diffusion evident in Figure 4.1c. 
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In this section the sensitivity of the computer model to 


changes in cell size (grid spacing) is investigated. The results 
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of a simulation using values of the parameters as listed in 
Table 4.1 are shown in Figures 4.2a to 4.2d. Figures 4.3a to 
4.3d are the results using parameters as in Table 4.1 except 
that INCX = 40 m and INCZ = 5m, i.e., the area of each cell 
was increased fivefold. 

Immediately obvious from a comparison of Figures 4.2a 
and 4.3a with 4.2b and 4.3b is the fact that the smaller grid 
produced a smoother particle distribution. Both cell sizes 
showed a tendency for accumulation of particles in the lowest 
portion of the valley in the region of reduced wind speed. 
Also, both showed tendencies for particle accumulation ina 
preferred ring near the axis of the plume. This appears to be 
the path taken by the majority of particles around the valley 
vortex and may be due in part to convergence which may be 
produced by the interaction of the advection and diffusion 
velocities but more likely to the fact that most particles are 
emitted at this height and not diffused substantially from it 
(see Appendix A). Also of interest are the observations that 
the particle distribution did not reach steady state after 2 h 
model time (as expected) and that, at 2 h, a number of particles 
had diffused to the opposite side of the valley. These Figures 
show another interesting effect of changes in cell size. Ina 
purely Eulerian mode, increased cell size generally ''smears out'' 
the concentration field, that is, distributes the pollutant 


over a wider area. However, the effect of increased cell size 


i ni beteth 2s eye bie al 79 saul : 
{ Send a 
o7 of.) 2s9uphF ee os 68.4 gonugiF wt 
7 2 a oe 
; josaxe 1.8 sldsl of 26 2 :stamsied pnizy ai fuden i 
[isa (46s To’ ore ad? , 9.) 0 2 = SIA) bee Ce 
« * 7 
ty i? basen Toar ¢ 
acta asso V7’ We noel Tet) eeoivde yiste anne i- 
| 
™ 5 j e 
biig tei tse srt? 76ers Joe 4) :j d&.3 bas doe aivlw abet 
3% my ntO8 ; Jj i ag (aoe 5 bsabong 
wee i * 
Paaweal 5 fold +. ois inwose Abt Yonehne) 6 beweta 
: : a 
bas ; shaq nt e/4 ni veliev 2 fefreg 
» on 
F un i iE ' 3 sbria3 Kepare A-od oataé 1% 
oo} 276 = hs m iy 2 Sc4 +89 0087 besieteig 
fiev sid teva 12 t GO, 6m sci yd neve) Aveg 
e 2 i t | i ve 
- ‘7 Ol 2367 nm 2 od Yen the xem 
ion i 
COlf2619In!) sad Ve ssaueed 
: 2&4 {Tot 2 vty i). ston dud Selidiggt 
: ’ : 
o7F Vitel ledud byzut fon ban tdeled eify te baeeiew 
i 
3 ? rs 7 os fh : Siar sei rh. s he) 
of Ci. - » " 7 a aw Th o hoe al 
. 1 4 ose ik ones ‘ a5 oie! tC slatiaeq eAD 
: : 
zaia} 2764. td° Jedh | ch @ <fady, Sap (be fosqne” 2a) sinkd aie. 
eOTURIT geet! i 2 BI) 20Qgo aniegs bezel tka 
6 ml «ete Ffss’ni @apients to Fostt= Bal- general ston woe < 
tt fu gt | » ~ ; j U ’ 
guo e169 Yi 2D Bxic. ites z62 130) ,sben nélial oa ott 
- wt! fT ve 1 ke : - ~ 
Itmtul lac on ead bI240 #1 3603 AY7 “fe ) S675 Sangd - on 


asic ites bazery 341s Io 328 
’ oS a | os : _ 
Sil. — ee iH a Seo 


t15, a3 eDavoue esa iy | 


50 


4074 


M) 


( 


307 


20-4 


planeta pil 


10-4 


0 100 200 300 400 500 600 700 800 
DISTANCE «ch 


Figure 4.2a: Particle distribution produced by standard parameters (Table 4.1) at 
t=3610s. 
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Figure 4.2b: Particle distribution: produced by standard parameters (Table 4.1) at 
t=718ls. 
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Figure 4.2d: Vertical profile of concentration 50m upslope from source 
produced by porameters.of Table 4.1. 
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Figure 4.3a: Particle distribution at t=3608s produced by parameters of Table 4.1, 
except that INCX=40m and INCZ=5m. 
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Figure 4.3b: Particle distribution at t=7195s produced by parameters of Table 4.1, 
except that INCX=40m and INCZ=5m. 
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Figure 4.3c: Vertical profile of concentration 100m downslope from source pro- 
duced by parameters of Table 4.1, except that INCX=40m and 


(Peo 


INCZ=5m. 
5 i= Sa 
5 Onno M 
uy a die SM Ly 
* 6.5 M 
1 3 
° 
Zz ry 
ean = © ame aes re 
eo 
© 
° 
° 
iu} 
° 
1 ° D L 4 
° 
® + a PS ere Cae ee ey: 
Pdi  -gunigus Shaae + : 
ey + 
r] + 
0 OE, oe Pe ee. a a 
0 600 1200 1800 2400 3000 3600 4200 4800 5400 6000 6600 7200 
WME USE! 


Figure 4.3d: Vertical profile of concentration 50m upslope from source produced by 
parameters of Table 4.1, except that INCX=40m and INCZ=5m. 
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on the present model appeared to be opposite. Particles 
accumulated near the upper boundary (and to a lesser extent 
near the valley center) rather than dispersed. The reason for 
this is unknown but indicated that the model was sensitive to 
changes in grid size. 

The vertical profiles of concentration at the downslope 
''station'' (100 m down the slope from the source), Figures 4.2c 
and 4.3c show interesting differences. While maximum 
concentrations in both cases were between 5 and 6 ppm, in the 
case of the larger grid the maximum occurred in the lowest cell, 
as might be expected. For the smaller grid the maximum occurred 
near 1.5 m. Again, this may be a result of the preferred track 
produced by the advection field, since movement by advection 
dominated movement by diffusion. Concentration curves were 
smoother for the larger grid because gain or loss of particles 
resulted in a smaller fractional change in concentration. 

Figure 4.2c illustrates that, for the smaller grid, concentrations 
in the lowest two cells were similar and variable for 
approximately the first hour. At that time the return flow had 
accumulated a sufficient number of particles in the cell 

centered at 1.5 m so that differentiation was possible 
Concentrations in the highest and lowest cells remained nearly 
equal and constant near 3 ppm. The concentration of the cell 
centered at 3.5 m in Figure 4.3c was much less than that of 

Figure 4.2c. Presumably this was due to its larger sampling 


interval in the vertical. The concentration of the 8.5 m cell 
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in Figure 4.3c was not influenced by the return flow. 

Figures 4.2d and 4.3d show the vertical profile of 
concentration 50 m upslope from the source. All effects here 
were due to return flow; diffusion from the (downwind) source 
did not occur. The maximum concentrations observed at the 
0.5 m and 3.5 m levels of the small grid case (Figure 4.2d) 
were approximately twice those of the large grid. The 
presumption that this is due to a smaller sampling extent in 
the vertical appears to be valid for the 3.5 m cell but not 
for the 0.5 m cell. Referring to Figure 4.2d, both the 0.5 m 
and 1.5 m cell reached concentrations of approximately 3.5 ppm. 
The 0.5 m cell of Figure 4.3d covers about the same vertical 
extent as the combined 0.5 m and 1.5 m cells of Figure 4.2d 
and yet the concentration was much less. The reason for this 
is not clear. 

Comparison of Figures 4.2d and 4.3d also shows that the 
onset of fumigation occurred slightly earlier in the small grid 
case. An explanation for this may be found in the influence 
of the grid size that appears in the calculation of the timestep 
by (3.13). A larger timestep, combined with poorer resolution 
of the velocity field, may account for the delay in the case 
of the larger grid. 

Another observation common to Figures 4.2d and 4.3d and 
reinforced by the results in Appendix A is that particles in 


the return flow arrived first in the upper-most of the lowest 
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three cells. This is consistent with the decreased path length 
and small vertical velocity gradient for elevated particles in 
the advection field. 

Near the valley slope the plume was typically 5 mor less 
thick. The vertical dimension of the cells chosen ranged from 
2 to 5m. Because of this, changes in grid size resulted in 
large changes in cell concentration. This implies that smaller 


vertical grid sizes are needed in this application. 


4.4 Sensitivity to Initial Source Size 


Recall from section 3.6 that the initial size and shape 
of the automobile exhaust plume was produced by a normal random 
number generator with given horizontal and vertical standard 
deviations. In this section sensitivity of the model to 
changes in the initial size of the exhaust plume is investigated. 
The basis of comparison is Figure 4.2, produced using the 
parameters of Table 4.1. Figure 4.4 was produced using the 
parameters of Table 4.1 except that SDEVX = 6 m and SDEVZ = 4 m. 
The initial size of the plume was doubled. The most important 
consequence appears to be the doubling of the initial vertical 
extent of the plume. 

Figures 4.2a, 4.2b, 4.4a, and 4.4b illustrate that, as 
expected, the plume produced with the larger initial size 


retained this size advantage throughout the 2 h simulation. 
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Figure 4.4a: Particle distribution at t=3610s produced by parameters of Table 4.1, 
except that SDEVX=6m and SDEVZ=4m. 
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Figure 4.4b: Particle distribution at t=7188s produced by parameters of Table 4.1, 
except that SDEVX=6m and SDEVZ=4m. 
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Figure 4.4c: Vertical profile of concentration 100m downslope from source pro- 
duced by parameters of Table 4.1, except that SDEVX=6m and 
SDEVZ=4m. 
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Figure 4.4d: Vertical profile of concentration 50m upslope from source produced 
by parameters of Table 4.1, except that SDEVX=6m and SDEVZ=4m. 
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Although the plume had a greater vertical extent in the large 
source case, the heights of the preferred path in both cases 
were similar. This again suggests that the preferred path was 
a result of the height at which the majority of particles were 
emitted. 

Figures 4.2d and 4.4d reveal that, as expected, the 3.5 m 
cell received a larger share of the particles in the case of 
the larger source size. These Figures show that the concentration 
in the 3.5 m upslope, large-source station reached more than 
2 ppm, double that of the small source case. There was a general 
shift to larger concentrations in the upper cells. In both 
cases maximum concentrations at the upslope position were near 
4 ppm at the end of the 2 h simulation. 

Similarly, a comparison of Figures 4.2c and 4.4c shows an 
increase in the concentration of the upper cell in the large- 
source case. The concentration increased from about 3.5 ppm 
to 4.5 ppm. The lower two cells both experienced decreases in 
concentration. In the 0.5 m cell, it dropped from approximately 
3 to less than 2 ppm. At the 1.5 m level, the concentration 
decreased to 4 ppm from about 5.5 ppm. Although the maximum 
concentrations were lower in the large-source case, the particles 
were distributed more evenly in the upper cells. This is 
evident in Figure 4.4b. 

One indication that the pollutant was more evenly distributed 
in the large-source case is shown by the uniformity of the 


1.5 m concentration. Figures 4.4c and 4.4d illustrate that, at 
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this height, concentrations at the upslope and downslope 
positions were approximately equal at most times. With increased 
source size, concentrations apparently became less dependent 


On position of measurement on the slope. 


4.5 Sensitivity to Particle Density 


This section examines the sensitivity of the model to a 
reduction in the particle density necessary for adequate 
resolution. Simulations were made using the parameters of 
Table 4.1, except that PWT = 0.4 g CO, i.e., compared to Figure 
4.2, approximately one-half the number of particles were created 
at each timestep, each representing twice the mass of CO. 
Figure 4.5 is the result of this simulation; Figure 4.2 is the 
basis for comparison. 

The particle position plots, Figures 4.5a and 4.5b, yield 
no new insights. There are no significant differences between 
these Figures and Figures 4.2a and 4.2b, except for the 
expected density difference. 

Comparison of Figures 4.2c and 4.5c reveals no change in 
average concentration in any of the cells. However, there is 
more variability in the curves representing the less dense 
particle distribution, Figure 4.5c. Similar comments are valid 


concerning comparisons of Figures 4.2d and 4.5d. 
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Figure 4.5a: Particle distribution at t=3609s produced by parameters of Table 4.1, 
except that PWT=0.4g CO. 
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Figure 4.5b: Particle distribution at t=7189s produced by parameters of Table 4.1, 
except that PWI=0.4g CO. 
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Figure 4.5d: Vertical profile of concentration 50m upslope from source produced 


by parameters of Table 4.1, except that PWT=0.4g CO. 
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It can be shown that, for the parameters of Table 4.1 at 
Standard temperature and pressure, each part per million of 
concentration corresponded to approximately 12 particles ina 
complete (not disected by the slope) cell. Therefore, each cell 
with a concentration of 3 ppm contained about 35 particles. 

In the case of the less dense distribution, 3 ppm resulted from 
about 17 particles. Thus, an increase in the variability of 
the concentration in the latter case might be expected but 

that adequate particles should exist to provide relatively 
smooth concentration curves. 

Although more smoothing can be applied to the curves 
output by the computer program, little can be done efficiently 
at intermediate timesteps. Because of the technique's hybrid 
quality, both the concentration field and particle distribution 
would need to be smoothed. Adjusting the latter is difficult 
and may rival a simple increase in density in terms of cost 


effectiveness. 


4.6 Sensitivity to Magnitude of Eddy Diffusivities 


4.6.1 No Diffusion 


This section compares the standard case (Figure 4.2) which 


used the parameters of Table 4.1 to the case for which no 


diffusion is permitted (advection only). The parameters for 
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this case are those of Table 4.1, except that DKMX =O and 
DKMZ = 0. Results are shown in Figure 4.6. 

Comparison of Figures 4.2a, 4.2b, 4.6a, and 4.6b reveals 
nothing startling. No particles appeared in the left-hand 
vortex of the valley in Figure 4.6b. Thus, this aspect of 
Figure 4.2b may be attributed correctly to diffusion. Particles 
in the diffusion case tended to greater accumulation near the 
bottom of the valley. The reason for this is not entirely clear. 
In the case of no diffusion, particles were advected parallel 
to the slope; very few particles were reflected by the slope. 
This was discovered during model development. When the diffusion 
process was included particles gained a velocity component 
normal to the slope. This could result in the following action. 
First, particles will tend to positions nearer the slope. Thus, 
they will be bound by pseudo-velocity streamlines which will 
carry them nearer the bottom of the valley (recall Figure 2.1). 
Advection velocities near the bottom of the valley are small. 
Therefore, interaction of advection and diffusion may allow 
particles to accumulate. 

Concentration profiles, Figures 4.2c, 4.2d, and 4.6c, and 
4.6d, at upslope and downslope positions reveal little of 
interest. Concentrations in the uppermost cell were generally 
slightly lower in the case of no diffusion - fewer particles 
reached that level without the aid of diffusion. At the 
downslope station the lack of diffusion appeared to result in 
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Figure 4.6a: Particle distribution at t=3604s produced by parameters of Table 4.1, 
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Figure 4.6b: Particle distribution at t=7181s produced by parameters of Table 4.1, 
except that no diffusion is allowed - advection only . 
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Figure 4.6d: Vertical profile of concentration 50m upslope from source produced 
by parameters of Table 4.1, except that no diffusion is allowed - 
advection only . 
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reasonable because, at the downslope station, the middle of 
the three cells is on the preferred route for particles 
circumnavigating the valley vortex (Figure 4.6b). Based on 
observations of concentration variability before averaging, 
the increase is probably not significant. Concentrations in 
the uppermost and lowest cell of Figure 4.6c appeared to be 
relatively uninfluenced by lack of diffusion. 

At the upslope position (Figure 4.6d), concentration in 
the lowest cell increased slightly in the absence of diffusion. 
Presumably, this is because the preferred route of the particles 
is nearer the slope at the upslope station. The differences 
in concentration with and without diffusion were small. The 
differences in height of maximum concentration between the 
upslope and downslope positions was probably an artificiality 
of the model advection field. The exclusion of diffusion from 
particle transport had little effect on concentrations within 


the parameter range investigated. 


4.6.2 A Range of Eddy Diffusivities 


In this section the sensitivity of the results to the 
magnitudes of the horizontal and vertical diffusivities is 
examined. No comparisons are made to Figure 4.2. Parameters 
are those of Table 4.1, except for cell size and variable 


diffusivities. In these comparisons INCX = 40 m and INCZ = 5 m. 
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In all cases horizontal diffusivity and vertical diffusivity 
were varied concurrently. Table 4.2 contains Figure numbers 


and corresponding values of eddy diffusivities for this section. 


Table 4.2: Eddy diffusivities and corresponding Figures 
for section 4.6.2. DKMX and DKMZ are, respectively, 
horizontal and vertical eddy diffusivity. 


Figure DKMX (m@s~!) DKMZ(m@s_ ") 
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As expected, Figures 4.3, 4.7 and 4.8 (a and b) show a 
general dilution in concentration throughout the valley vortex 
as diffusivities increase. In addition, as diffusivities 
increased more and more particles were diffused into the 
adjacent vortex. This accounts for part of the dilution. By 
inspection only, it is difficult to see any increase in plume 
spread caused by increasing diffusivities. Again, the diffusion 
process was nearly masked by advection with diffusivities of 


the magnitudes used here. 


Figures 4.3c, 4.7c and 4.8c compare the results at the 


downslope station. Concentrations in the 8.5 m and 3.5 m 
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Figure 4.7a: Particle distribution at t=3617s produced by parameters of Table 4.1, 
except that INCX=40m, INCZ=5m, DKMX =1.0x 10°'m’s',and 
DKMZ=5.0 x 10‘m’s'. 
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Figure 4.7b: Particle distribution at t=7204s produced by parameters of Table 4.1, 
except that INCX=40m, INCZ=5m, DKMX= 1.0x 10 m's',and 


DKMZ=5.0 x10 ms 
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Figure 4.7c: Vertical profile of concentration 100m downslope from source pro - 


duced by parameters of Table 4.1, except that INCX= 40m, INCZ=5m, 
DKMX=1.0x 10 ms. rand DKMZ=5.0 x10 m’s*. 
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Figure 4.7d: Vertical profile of concentration 50m upslope from source produced 


by parameters of Table 4.1, except that INCX=40m, INCZ=5m, 
DKMX=1.0x 10°m's',and DKMZ=5.0x 10 ms. 
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Figure 4.8a: Particle distribution at t=3601s produced by parameters of Table 4.1, 
except that INCX=40m, INCZ=5m, DKMX=1.0x 10 ms .and 
DKMZ=5.0 x 10 m’s' 
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Figure 4.8b: Particle distribution at t=7218s produced by parameters of Table 4.] 
-| é 
except that INCX= 40m, INCZ=5m, DKMX=1.0x 10 m’s ,and 
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Figure 4.8c: Vertical profile of concentration 100m downslope from source pro- 
duced by parameters of Table 4.1, except that INCX=40m, INCZ=5m, 
DKMX=1.0x 10'm’s',and DKMZ=5.0 x10 m’s". 
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Figure 4.8d: Vertical profile of concentration 50m upslope from source produced 
by parameters of Table 4.1, except that INCX=40m, INCZ=5m, 
DKMX=1.0x 10 m s',and DKMZ=5.0 x10 ms. 
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cells were relatively uninfluenced by changes in diffusivities. 
Concentrations in the lowest cell decreased with increasing 
diffusivities. With increased diffusivities, particles were 
diffused away from the valley slope. This was the expected 
result. The changes in concentration resulting from a 100-fold 
change in magnitude of the eddy diffusivities were relatively 
small, less than 1 ppm (about 20%). Maximum concentrations, 
found in the lowest cell, ranged from about 4 to 5 ppm. 

Concentrations in the lowest cell at the upslope station 
(Figures 4.3d, 4.7d and 4.8d) show much the same variation 
with diffusivities as those observed at the downslope station. 
Concentrations in the lowest cell varied by approximately 202, 
ranging about the value 2 ppm. Unlike the downslope case, 
concentrations in upper cells also varied somewhat with 
diffusivity. In these cells concentration increased with 
increasing diffusivities. Concentrations in the 8.5 m cell of 
Figure 4.8d show a large increase over those in Figure 4.3d. 
This may indicate that effective velocities produced by 
diffusion rivalled those of advection at this height. 

Results of this section, and those of section 4.6.1, 
suggest that concentrations are relatively insensitive to 
large changes in the magnitude of both horizontal and vertical 
eddy diffusivities. Magnitudes of the diffusivities (if not 
their spatial and temporal variation) investigated here seem 
reasonable for the strong inversion and large horizontal 


temperature gradient observed in the North Saskatchewan River 
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valley in Edmonton. Because concentrations were relatively 
insensitive to changes in diffusivities, this suggests that 
the valley advection wind regime plays an important, if not 


dominant, role in determining pollutant levels in the valley. 


4.7 Sensitivity to Advection Velocity 


This section examines the sensitivity of the model to a 
reduction by a factor of two of the magnitude of the advection 
velocity. Parameters used are those of Table 4.1, except for 
cell size and magnitude of advection velocity. For this 
analysis INCX = 40 m, INCZ = 5 m, and AVWF = 0.25 aa. Figure 
4.3 is the basis of comparison. Figure 4.9 shows the results 
for a reduced wind speed. Figures 4.9a and 4.9b reveal that, 
as expected, the distance travelled by the particle plume was 
much smaller when the wind speed is reduced. In addition, the 
plume appears more dense. The reduction of wind speed has 
constrained the same total number of particles to a smaller 
volume. 

Figures 4.3c and 4.9c indicate a substantial increase in 
concentration in the lower levels at the downslope station in 
the case of reduced advection velocity. Return flow of 
pollutants did not dramatically increase the concentration in 
the 0.5 m cell, at least in the time period modelled, in the 


case of reduced wind speed. Instead the concentration was 
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Figure 4.9a: Particle distribution at t=3633s produced by parameters of Table 4.1, 
except that INCX=40m, INCZ=5m, and AVWF=0.25ms’. 
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Figure 4.9b: Particle distribution at t=7225s produced by parameters of Table 4.1, 
except that INCX=40m, INCZ=5m, and AVWF=0.25ms.. 
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Figure 4.9 c: Vertical profile af concentration100m downslope from source produced 
by parameters of Table 4.1, except that INCX=40m, INCZ =5m, and 
AVWF=0.25ms_. Concentrations at 0.5m for the period 900s to 1800s 
are greater than 7.5 ppm. 
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Figure 4.9d: Vertical profile of concentration 50m upslope from source produced 
by parameters of Table 4.1, except that INCX=40m, INCZ=5m, 


and AVWF=0.25ms: 
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maintained at a high level. The 0.5 m curve of Figure 4.9c 
illustrates the decrease in source strength, at least during 
the first hour of model time. As expected, the increase in 
concentration with reduced advection velocities was less at 
greater heights. The 8.5 m cell showed negligible change in 
concentration. 

Figure 4.9d shows that two hours was insufficient for 
investigation of concentration changes at the upslope position 
when the slope wind was reduced. The concentration in the 
uppermost cell was substantially larger in Figure 4.9d than in 
Figure 4.3d. This may be because the plume had not progressed 
as far into the upper right corner of the valley (Figure 4.9b). 
The plume path was somewhat foreshortened; the descending region 
of the plume was nearer the upslope station in the case of 
reduced wind speed. Concentrations in the lower cells increased 
rapidly at the end of the simulation. The fastest particles 
traversed the vortex in about 3600 s, double that of Figure 4.3d. 
This was as expected. 

The results of this section suggest that concentrations 
in the valley are sensitive to changes in the magnitude of the 
advection velocity. This was an expected result. However, 
observations in the North Saskatchewan River valley (Hage, 1979) 
indicated that, once developed, the slope wind was steady until 
its demise at sunrise. Thus, the results of this sensitivity 
to wind speed analysis may be of little use in elucidating 


aspects of real valley winds but do indicate that the model 
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predicts the correct trend of increasing concentrations with 


decreasing wind speed. 


4.8 Computer Requirements 


The computer code of this model was written for and run 
On the University of Alberta's Andahl 470V/6 computer. The 
program was written in FORTRAN IV and compiled using the 
FORTRAN H compiler. The University's Andahl computer at 
present has an 8 M byte core storage capacity. The Michigan 
Terminal System (MTS) operating system used at the University 
of Alberta provides large amounts of virtual memory. 

Central processing unit (CPU) time and storage space 
required by the code were dependent on the application. The 
standard run, using the parameters of Table 4.1, required 
approximately 100 s CPU time. Reducing by a factor of two 
the number of particles produced at each timestep (section 4.5) 
reduced CPU time by about 40% to about 60 s. Reducing the 
area of each cell by a factor of five (section 4.3) reduced 
CPU time by 60% to approximately 40 s. Reducing the maximum 
velocity on the grid by a factor of two also halved CPU time. 
Changing other parameters had little effect on CPU time. Time 
required was found, at least in the cases examined, to be 
primarily dependent on three parameters: vertical grid spacing; 


vertical velocity and; vertical eddy diffusivity. Combinations 


=| 


4 


Oi SMES oo. ae3 7 
bon laxe teea 
spNnisuqe biND heat dyay 


enory 1) 1 amso sya iviau 


y ’ 


e'erisdtA/ te ysierawint ond ne 
Vi VARTSOF Gf cardtinw see rare 
fesavthw anys 400% a3 AEA 
"tle atts <¢0c 4 8 ons 26d eee ze 

ot iaado: (2TH) we eve Isaianets 
‘ITNUONe so rst 2sbhivosg coved tA. es - 
4) ) oni 2zenorwg ‘eee @ : 


ry 1% =a" sins ade 7d ug toper : 
vn 


Né1609 SAT onieb , nia OI SDOGPE! 7 


Ps én » ¥ 4 a 
. i(9 299 »« DG VisTenteargqee 


. 


se’ tw 2 4 wr) ome hs a" - 
as ‘ € _— 7 USg teswiot 7 


- ~ 
Oo cLiotG 2elo isco Ghee aes 
7 ‘ 


c toi3s* 5 vd Miso ieee ae ae 
pn 
‘ : g } 
Y' 5 e*t no ges ' oF Je on wnt s 


is ya tisp tt wo letige 


. c 


_ 


ine 
aa 


ii! bet eqetiire tag orge! wre 


At trasl-gy «Dag #6 bavh a 
16159 sei) aa nel oe “ aba = 
fox 


* (noite Bie vis 


any 


87 
of decreased vertical grid spacing, increased vertical velocity, 
Or increased vertical eddy diffusivity tended to decrease the 
length of the timestep and so increase CPU time. 

Storage space required by the compiled version of the 
computer program ranged from 0.4 M to 0.5 M bytes, dependent 
mostly on cell size and number of particles created each timestep. 
In all cases, the charge for CPU time was approximately equal 
to the charge for storage space. Total cost for each run (not 


including plotting) was approximately double the cost of the 


CPU time alone. 
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CHAPTER V 


SUMMARY AND CONCLUSIONS 


5.1 Summary 


A series of micrometeorological experiments were conducted 
in the North Saskatchewan River valley in Edmonton during 1977 
and 1978. These experiments provided data from which the 
microclimatology of small valleys may be deduced. Resultant 
information about the occurrence of pollution episodes and the 
levels of pollutant during these episodes may be of interest 
to those who plan the use of city valleys. This thesis is one 
part of the effort to understand the micrometeorology of smal] 
urban valleys. It is a first step in modelling the transport 
of pollutants within these valleys. 

A two-dimensional model was formulated using a particle- 
in-cell approach. Particles, representing specified amounts 
of CO, were created within a V-shaped model valley. These 
particles were advected by a wind field that was specified for 


the model, and dispersed by a simulated turbulent diffusion 
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process. The model used a concentration-gradient diffusion 
technique for distributions of particles large compared to the 
grid size. For smaller distributions, diffusion was assumed 
to be Gaussian in character. The two were found to produce 
approximately equal rates of diffusion. The source of the 
particles was a point on the valley slope, representing an 
infinite line source which in turn represented a roadway 
running along the valley axis. The advection field was a 
(somewhat artificial) double vortex representing potential flow. 
Particles were allowed to diffuse, but were not allowed to be 
advected, to the adjacent vortex and out of the top of the 
valley. No interaction of the valley pollutant distribution 
or wind field with those of the overlying flow was allowed. 
The valley was essentially a closed system. 

The model depicted a situation in which an inversion, and 
thus the double vortex valley wind system, filled the entire 
depth of the valley. The slope wind speeds of the advection 
field were similar to those that have been observed in the 
North Saskatchewan River valley in Edmonton. Eddy diffusivities 
used by the model were Fickian, i.e., constant in space and 
time. Arguments for approximating the magnitudes of the 
diffusivities were based on the observed intensity of inversions 
and horizontal temperature gradients deduced to exist near the 
slope. 

The model correctly predicted an increase in concentration 


with decreasing wind speed and a decrease in concentration with 
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decreasing source strength, although these results were masked 
to some degree by the effects of the closed valley system. 

These effects were responsible for the model concentrations 

not reaching steady state. Concentrations downslope from the 
source reached over 6 ppm. The return flow of the double vortex 
advection field was capable of producing concentrations of up 
to 5 ppm. This concentration was found somewhat upslope from 
the source, out of the range of direct diffusion. 

At the upslope position, the effects of the return flow 
were first observed in cells somewhat above the slope. Some 
minutes later, concentrations at the slope also increased. 

The fastest particles traversed the vortex in approximately 

25 minutes in a flow with downslope wind speed of 0.5 oe 

For a wind speed of approximately 0.25 ae travel time 
increased to about | hour. In most sensitivity analyses, 
concentrations at the downslope station appeared to reach their 
maximum values within the two hour simulation. Concentrations 
at the upslope station appeared to be increasing at the end of 
the simulation. Concentrations near the slope were found to 

be relatively insensitive to the magnitude of eddy diffusivities 
and to changes in particle density. Increasing the initial 
size of the source generally decreased concentrations downslope 
from the source. Concentrations at the upslope station were 
unchanged by changes in initial source size. At both stations 
increasing the initial source size decreased concentrations at 


the slope and increased concentrations above the slope, but 
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with no changes in maximum observed concentration. Decreasing 
by a factor of two the magnitude of the advection velocity 
doubled the concentration downslope from the source during the 
first hour and increased it by about 20% at later times. 
Decreasing the advection velocity by a factor of two doubled 
the transit time of particles. Two hours was an insufficient 
time to investigate concentrations at the upslope station in 
this case. "increasing by’a factor ‘of five ‘the arearof “each 
cell resulted in generally lower concentrations at similar 
heights above the slope. Somewhat surprisingly, the particle 
scatter graph for this case suggested an increase in 
concentration in some parts of the valley, at odds with the 
usual result found when using Eulerian methods. The reason for 
this was not determined but it is suggested that cell sizes 
were too large. 

It was found that magnitudes of vertical components of 
model parameters (i.e. cell size, wind speed, and eddy 
diffusivities) were the limiting factors in determining the 
cost of computation. This was due principally to the scale of 


the model valley. 


Bi Zeeconorusions 


Within the limits of the present study (small valley, 


closed valley system, and single, along-valley line source) 
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the ground-level concentrations downslope from the source were 
controlled by the interaction between advection and initial 
source diffusion. Concentrations upslope from the source 
apparently were influenced less by initial source diffusion. 
Atmospheric diffusion following the initial ne growth appeared 
to be of secondary importance for diffusivities that were 
believed to be realistic in moderate to intense valley inversions. 
The time required for a particle to traverse one circuit 
of the valley vortex appeared to be determined predominantly 
by advection when diffusion was represented by realistic 
diffusivities. Although cycle times were found solely from 
the time of first arrival at the upslope cell, they were similar 
to the cycle times suggested by Paterson and Hage (1979) when 
the speed of the downslope wind was about 0.25 cell Cycle 
times of one hour and typical downvalley wind speeds of 0.5 7s 
suggested the along-valley length of one coil of the helix to be 
about 2 km. Therefore, in down-valley winds of this magnitude 
or greater, recirculation of the air past an along-valley line 
source will be important only if that source extends several 
kilometers upvalley. 
Results of investigations of the effects of changes in 
cell size on predicted concentrations suggested that the cel] 
sizes used in this application were too large. This throws 
some doubt on the validity of the results. Since the vertical 


dimension of the plume was about 5 m, vertical grid size 
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should be substantially smaller, say 0.5 m. Because vertical 
components were the limiting factor, computer costs would rise 
accordingly. 

Results showed that predicted concentrations were 
sensitive to initial source diffusion. This lends Support to 
the conclusion that diffusion following the initial plume 
growth is of secondary importance. It also Suggests that more 
accurate parameterization of the source is needed. In particular 
this entails a better representation of the effects of buoyancy 
and mechanical mixing and perhaps a more detailed analysis of 
the vehicle emission factor. 

The magnitude of the concentration in the return flow is 
important to discussions of the possible existence and structure 
of such a flow in small valleys. A typical concentration near 
the slope surface upslope from the source was 4 ppm after two 
hours. During this time source strength decreased by nearly 
40%. Concentrations of 4 ppm support the existence of return 
flow. The value is substantial; model refinements such as 
interaction with the overlying flow and possibly larger initial 
plume growth would serve to decrease the value but presumably 
not to an undetectable size. The concentration profiles do 
little to indicate the depth of the slope wind. The depth of 
the plume resulting from diffusivities which are believed to 
be realistic indicates, however, that flow in the lowest few 


meters is important for determining concentrations near the 


source. 
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In conclusion, the PIC technique appears to be of value 


in modelling pollutant transport in small valleys with return 


flow. 


5.3 Suggestions for Future Study 


The basis of the model - representing vehicular emissions 
by Lagrangian particles in an Eulerian grid - appears to have 
potential for modelling transport in small valleys. However, 
it is important to develop more realistic boundary layer and 
source parameters. Although the model wind field appears to 
simulate adequately the slope wind, it is thought to be less 
realistic at the top of the valley. The addition of an 
overlying potential flow with appropriate interaction at the 
interface is a possible solution. This should allow for some 
pollutant dilution. 

Parameterization of the eddy diffusivities needs further 
refinement. In the moderate to intense inversions found in 
small valleys, the effects of mechanical mixing become 
important. The region beneath the level of the maximum slope 
wind could be characterized as of neutral stability, as 
suggested by di Cenzo (1979), thus increasing dispersion in 
the lowest levels. 

Initial source diffusion requires more realistic 


parameterization. In this application buoyancy was ignored. 
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Because of the subordinant role of diffusion at locations away 
from the source, a change in the effective height of the plume 
by only a few metres may be important for predicting 
concentrations. If, as suggested, mechanical mixing is 
important, then perhaps the initial size of the plume also needs 
to be increased. 

The source inventory could be suitably enhanced. Drainage 
of pollutants from the city might be simulated by positioning 
sources of appropriate strength at the upper slopes. Upvalley 
sources could be accounted for by assuming a background additive 
concentration. 

Finally, the sensitivity of the model to cell size should 
be investigated further. In particular, cells that are smaller 
than those employed here should be used. In addition, the 
trajectory of several representative particles should be 
determined. This may clarify typical traverse times for particles 


in the valley vortex. 
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APPENDIX A 


AN ILLUSTRATION OF A PARTICLE SCATTER 


DIAGRAM DISPLAY 


This appendix contains a series of particle scatter graphs 
at intervals of 900 s. The parameters used in these simulations 


Be and 


are those of Table 4.1, except that DKMX = 1.0 x 10. 
DKMZ = 5.0 x 10 -més |. 
The dashed line in Figures Al to A5 is an advection 
velocity streamline. It is meant to pass through the center 
of mass of the source puff. For a continuous half-plane Gaussian 
distribution, the center of mass is near 0.680. For a puff 
with o = 2 m, the center of mass should be near 1.4 m. Thus, 
the streamline in the Figures is near the center of the initial 
source configuration. At early times (Figures Al and A2), it 
is obvious that those particles created in the upper levels 
of the plume traverse the valley vortex faster than those 
created in the lower levels, due to shorter path length. 
Particles require a long period of time to reach the bottom 
of the valley because of the rapidly decreasing velocities 
near the valley bottom. 


It might be expected that, at later times, the illustrated 


streamline remains near the center of the plume. Figures A3 
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to A5 show that this is not the case. Deviations of the 
streamline from the center of the plume may be caused by 
the presence of diffusion or by errors produced by inadequate 


resolution in the finite-difference representations. 
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COMPUTER PROGRAM 


This appendix contains the computer code on which this 


study is based. 
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WRITTEN BY: R.RUDOLPH 


PURPOSE: TO MODEL THE TRANSPORT OF POLLUTANTS IN A VALLEY 


THIS PROGRAM MODELS THE ADVECTION AND DIFFUSION OF AEROSOLS AND IS 
BASED ON THE WORK OF ROLF LANGE (REF: JAM( 1978), 3, PP.320-329).THE 
REFERENCE CONTAINS A MORE COMPLETE BIBLIOGRAPHY. 


THE CODE SOLVES THE TWO DIMENSIONAL ADVECTION - DIFFUSION EQUATION BY 
THE PSEUDO - VELOCITY TECHNIQUE FOR A GIVEN NON-DIVERGENT ADVECTION 
FIELD. THE METHOD IS BASED ON THE PARTICLE - IN - CELL TECHNIQUE WITH 
THE POLLUTANT CONCENTRATION REPRESENTED BY LAGRANGIAN PARTICLES IN AN 
EULERIAN GRID MESH. 


THE FOLLOWING FORTRAN UNITS ARE USED: 


UNIT 3 : CONCENTRATION OF ONE CELL AS A FUNCTION OF TIME 
UNIT 4 : STREAM FUNCTION ARRAY 

UNIT 5 : INPUT DATA (SEE BELOW) 

UNIT 6 : DIAGNOSTICS 

UNIT 7 : DOCUMENTATION OF ERRORS 

UNIT 8 : SOURCE UNIT (*SOURCE*) 

UNIT 12 : PARTICLE POSITIONS 

UNIT 13 : CONCENTRATION ARRAY 


THE FOLLOWING SUBROUTINES ARE CALLED FROM THE MAIN PROGRAM: 


SIRIMMEIN 2 (EANECQUIL/NMES, TWirle INE Ole (a/Nelri CEILI. rll (\ (leit alolipNILILAY 
BELOW GROUNDLINE,THE AREA IS SET TO -1. THE BOUNDARY 
TS @ Asc UMEDET CRE Gar Am Viseoeesit APE Dav AIaE Esvan 
ALTHOUGH MORE COMPLICATED VALLEY SHAPES MAY BE USED, 
IT IS ASSUMED THAT THE SLOPE IS CONSTANT ACROSS ANY 
ISEIL 


ADVECT : COMPUTES A STATIONARY ADVECTION FIELD DEFINED AT 
CELL CORNERS. 


SOURCE : DETERMINES THE NUMBER OF PARTICLES AND THEIR INITIAL 
POSITIONS PRODUCED AT EACH SOURCE LOCATION AT EACH 
TIMESTEP. THE DISTRIBUTION IS ASSUMED TO BE APPROX- 
IMATELY GAUSSIAN. 


SUM : DETERMINES THE PARTICLE CONCENTRATION IN EACH CELL 


DIEAe SAGAECULATESS A) DIFFUSTONS VEEOC iY SAT EACH CELE CORNER 
FROM THE CONCENTRATION GRADIENT. VELOCITIES ARE SET 
TO ZERO FOR EACH CORNER BELOW GROUNDLINE. 


INTER : INTERPOLATES ADVECTION AND DIFFUSION VELOCITIES 
FROM CELL CORNERS TO PARTICLE POSITIONS. FOR SUB 
GRIDBSCABEMoEZEDEPURRc OU pEUSTONSVELOGI IE Sm AT) 
PARTICLE POSITIONS ARE COMPUTED DIRECTLY FROM 
GAUSSIAN DISPERSION. THE TRANSPORT DISTANCE IS 
FOUND AND THE LENGTH OF THE NEXT TIMESTEP IS CALCU- 
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EATREDS 

TSTEP : ADVANCES PARTICLES ALONG PSEUDO VELOCITY STREAMLINES. 
IF NECESSARY, PARTICLES ARE REFLECTED AT THE GROUND- 
LINE. THE VARIANCE OF EACH PUFF IS CALCULATED. 


RENUM : PARTICLES TRANSPORTED OUTSIDE THE GRID ARE DESTROYED; 
THOSE REMAINING ARE RENUMBERED. 


OUTPUT : WRITES PARTICLE POSITION AND CONCENTRATION DATA 
INTO FILES APPROPRIATE FOR USE IN PLOTTING ROUTINES 


THE FOLLOWING SYSTEM/IMSL SUBROUTINES ARE CALLED: 
TIME : STORES THE NUMBER OF SECONDS SINCE 1900 A.D. IN ‘ITM’. 


GGGNOQF : DETERMINES ONE NORMALLY DISTRIBUTED RANDOM NUMBER. 


DEFINITION OF PRIMARY ARRAYS AND VARIABLES: 


rte s iinlz ARIE (iP le/\elrl @lsilit 

CONC : THE NUMBER OF PARTICLES/UNIT AREA OF EACH CELL 
UDR EC ROSSEVABEEN sO pew S LONMVIEIO Cie in, 

VDP ae VERT CALSDIEBUSTON VELOCITY 

VAR eECROSSe VALE EV ADV EGREONSVELOGC Ih. 

VA :; VERTICAL ADVECTION VELOCITY 

XP : CROSS VALLEY (HORIZONTAL) PARTICLE COORDINATE 

Ze VeRLCAME CART GUE GOORDINATE 

JFLAG : POSITION ARRAY. ie eLiVesi., IAN WEES ile) lls ile silos) 
DEEX™: = GROSS, VAEEEY TRANSPORT DISMANGE 

DEEZ VER TLCA RST RANSPORT SDS TANGE 

VAR : VARIANCE OF EACH PUFF 

NIT : NUMBER OF PARTICLES IN EACH PUFF 

NPIT : NUMBER OF PARTICLES PRODUCED PRIOR TO EACH PUFF 
Cie ae CREAT LONE TIME SO RSEACH SURE 

CMX : HORIZONTAL CENTER OF MASS COORDINATE OF EACH PUFF 
CMZ : VERTICAL CENTER OF MASS COORDINATE 

XS : CROSS VALLEY COORDINATE OF EACH POINT SOURCE 

INCX : SIZE OF GRID IN HORIZONTAL 

WING = WiaisrineiVk iia) Sil7dle 

HPRIME : DEPTH OF VALLEY 

H : UPPERMOST EXTENT OF THE GRID 

i 3 MAD) Ole Wale WANLILIEN 

NP : RUNNING TOTAL OF PARTICLES PRODUCED 

USTAR : FRICTION VELOCITY (INPUT) 

STAB : DIMENSIONLESS HEIGHT (Z/L); AN INDICATOR OF STABILITY 
NT : THE NUMBER OF PUFFS PRODUCED 

DELT : THE TIMESTEP INCREMENT 

1 2 WON, Witte (SIL/NPSIEls ihe SWhi le ALE Ile’ Ss 

NS : THE NUMBER OF POINT SOURCES IN THE GRID 


LOGICAE* 1 LEMT(1)/7 *7 / 

LOGICAL*4 INT ,MASK/ZOOOOF FFF / 

REAE* Seo E ED 

EQUIVALENCE (INT,ITM) 

DIMENSION AREA(41,26),CONC(41,26),UD(41,26),VD(41,26), 
1UA(41,26),VA(41, 26) 

COMMON XP( 10000), ZP( 10000) , JFLAG( 10000), 

1{DELX( 10000) , DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(SO00) ,CT(5000), 
4CMX (5000) ,CMZ(5000) ,XS(10), INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 


TOE sli Nis 
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INPUT VALLEY AND GRID INFORMATION: 


L : VALLEY WIDTH (IN METERS) 

H : UPPER BOUND (IN METERS) 

HPRIME : HEIGHT OF VALLEY RIDGE (BOTTOM AT (L/2,0)) 
INCX,INCZ : HORIZONTAL AND VERTICAL GRID SIZE (IN METERS) 


READ(5,LFMT) L,H,HPRIME, INCX,INCZ 


INPUT TIME INFORMATION 


DELT : LENGTH OF INITIAL TIMESTEP (LATER STEPS CALCULATED) 
TMAX : MAXIMUM TIME MODEL IS TO RUN (IN SECONDS) 

NTMAX : MAXIMUM NUMBER OF TIMESTEPS ALLOWED 

NS : NUMBER OF POINT SOURCES 


READ(S,LFMT) DELT,TMAX,NTMAX,NS 
INPUT SOURCE AND WIND INFORMATION: 


SDEVX,SDEVZ : INITIAL STANDARD DEVIATION OF EXHAUST PLUME 

WEB 2 WMS (EID DIFFS UNIOSS CONSIWANNR Ilr WES, weno 

DKMX ,DKMZ : VALUES OF CONSTANT DIFFUSIVITIES 

TAV : ARE ADVECTION VELOCITIES ZERO (DIFUSION ONLY) ? IF YES, IAV=0 
AVWF : ADVECTION VELOCITY WEIGHTING FACTOR IF AVWF=5, U=0(.5M/S) 
USTAR,STAB : FRICTION VELOCITY (M/S) AND STABILITY (Z/L) FROM DATA 
XS(NS) : X-POSITION OF SOURCES (IN METERS) 


READ(5,LFMT) SDEVX,SDEVZ,ICD,DKMX,DKMZ,IAV,AVWF ,USTAR, STAB, (XS(I) 
1,1=1,NS) 


INPUT OUTPUT INFORMATION: 


IPSO : IF IPSO=O, VALUES OF STREAM FUNCTIONS ARE WRITTEN ON UNIT 4 

NO : WRITE CONC. AND PART. POSITIONS ONTO UNITS 12 AND 13 EVERY NO SEC. 
CC) 2 WROTE CONC. Ol ONE SIPEG IIED GEIL (AVERY Wee) IWietESvleleS ONive) Whey 3s 
X,Z : X AND Z COORDINATES (IN METERS) OF ONE SPECIFIED CELL 


READ(5,LFMT) IPSO,NO,ICO,X,Z 


DATA NOI/O/,NOTIN/O/,1X/41/,12/26/ 
CNL Wuhe( WS .@, Wu) 
INT=INT.AND.MASK 
SEED=ITM 
CALL SRAREA(IX,1Z,AREA) 
CALL ADVECT(IAV,AVWF,IPSO,IX,1Z,UA,VA) 
© 
NT=O 
NP=O 
1OO3T=h+DEET 
CALL SOURCE(SDEVX,SDEVZ,SEED,IX,1Z) 
CALL SUM(ICO,X,Z,IX,1Z,CONC,AREA) 
CALL DIF1(ICD,DKMX,DKMZ,IX,1IZ,UD,VD,CONC) 
GALE INTER(SDEVX, SDEVZ,1CD,DKMX ,OKMZ, IOUT, IX,IZ,UD,VD,UA,VA) 
CALL TSTEP(IOUT,NOTIN) 
CALL RENUM 
eC. La Now) <e@ vo 2oao 
CALL OUTPT(IX,1IZ,CONC) 
NOI=NOI+NO 
200 IF(T.LT.TMAX.AND.NT.LT.NTMAX) GO TO 100 
CALL OUTPT(1X,1Z,CONC) 
STOP 
END 
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SUBROUTINE SRAREA(IXX,12Z,AREA) 


WRITTEN BY:R.RUDOLPH 


COMPUTES CELL AREAS FOR A GIVEN VALLEY SHAPE 


{ 


DIMENSION AREA(IXX,12Z) 


COMMON XP( 10000) ,ZP( 10000) , UFLAG( 10000), 
1DELX( 10000) ,DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(5000) ,CT(5000). 
1CMX (5000) ,CMZ(5000) ,XS( 10), INCX, INCZ,HPRIME,H,L,NP,USTAR,STAB.NT, 


DEET ENS 
PI=3.1415926 
S=2*HPRIME/L 


IGRIDX,IGRIDZ ARE NUMBER OF CELLS IN HORIZONTAL, VERTICAL 


IGRIDX=L/INCX 
IGRIDZ=H/INCZ 


DETERMINE CELL CORNERS 


IFLAG=0O 


DO AO w=). veR wz 
Z1=(1Z-1)*INCZ 
Z2=Z1+INCZ 


DO 190 IX=1,IGRIDX 
Med = Gl Xe) Ss NCX 


X2=X1+INCX 


N=0 


IF(Z1.GE.HPRIME) GO TO 50 


CHECK WHICH CELL CORNERS ARE ABOVE THE GROUNDLINE 


1 


A1=L/2-Z1/S 

POX) Gas lL/2) 
A2=L/2-Z2/S 

IF (X1.GE.L/2) 
A3=(L/2-X1)*S 
F(OC) Ele (L772) 
A4=(L/2-X2)*S 
THO GEeIl/ 2) 


At=L-A1 


A2=L-A2 


A4=-A4 
LE CC CL/29GE 21) FAND SOX: 


AND 1 (X4207.A1))) N=N+4 


Tia Ce 2G Ee) ANDi (Oxo 
PANDNGXOe ADD aN=N+ 10 
ite (CU £2 GE . Xe) AN. (Oe. 


.AND.(X2.LT.A2))) N=N+100 


TiBGGCEA2eGE axa) SANDE Oxaie 


.AND.(X1.LT.A2))) N=N+1000 


CHECK N AND CALCULATE APPROPRIATE 


20 
30 


IF(N. 
IF(N. 
IF(N. 
IF(N. 
IF(N. 
IF(N. 
IF(N. 
IF(N. 
IF(N. 
WRITE(7, 
FORMAT(5X,’SHAPE NOT ALLOWED FOR CELL’ ,2X,13,’,’,13,5X, ‘N= 


EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQR 
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1000 ) 
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60 
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110 
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AREA(IX,1Z)=-2 
CORTON 1S 

40 AREA(IX,1Z)=-1 
GORTOMSoO 

50 AREA(IX,1Z)=INCZ*INCX 
CORTOM IS © 

60 AREA(IX,1Z)=INCZ*INCX-((A1-X1)*(A3-Z1)/2) 
GOMICHA SO 

70 AREA(IX,1Z)=INCX*INCZ-((X2-A1)*(A4-Z1)/2) 
GORTOM1S© 

80 ECA LT Aa) EO TO Go 
AS=A3 
DIF=A3-A4 
(e{8) 1/8) KC) 

30 AS5S=A4 
DIF=A4-A3 

100 AREA(IX,12Z)=INCX*(Z2-A5)+INCX/2*(DIF) 
GO) TO) 150 

110 Z1INT=A1 
Z2INT=A2 
RC ZUNE [Oy ZOU) Co We) 2e 
AREA(IX,1Z)=(X2-Z41INT)*INCZ+INCZ*(Z1INT-Z2INT)/2 
GORTO mS © 

120 Z1INT=A1 
Z2INT=A2 
We (Vu Cl. ZANT) CO TO Ae 
AREA(IX,1Z)=INCZ*(Z1INT-X1)+INCZ*(Z2INT-Z41INT)/2 


GOMmIOMIS© 

130 AREA(IX,1Z)=(A2-X1)*(Z2-A3)/2 
GO TO 150 

140 AREA(1X,1Z)=(X2-A2)*(Z2-A4)/2 


150 RECT=INCX*INCZ 
IF(RECT.LT.AREA(1IX,1Z).OR.AREA(1IX,1IZ).EQ@.0) GO TO 160 
IF(AREA(IX,1Z).EQ.-1.0R.AREA(1IX,1IZ).GT.O)GO TO 190 
IF(AREA(IX,1Z).EQ.-2) GO TO 180 


160 WRITE(7,170)IX,1Z,AREA(IX,1Z),N 

ZO roledmvelap (She, (NRE le lEILIL? te, 957 USO. “Woy UN) PikoPlale Ry\Nele 
1a ON, TINPRIENG ? 1}. 2M, ONE? , We) 

180 RAG = A 

ie) CONTINUE 


200 CONTINUE 
i Ch RlaAGeNEh) mG Oi OMm220 
WRITE(8,210) 
210 FORMAT(//3X,’PREMATURE TERMINATION IN SRAREA. CHECK ERROR FILE.’ 
/ fa) 
STOP 
220 RETURN 
END 
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SUBROUTINE ADVECT(K,C,KKK,1IXX,1ZZ,UA,VA) 


PENS EY aoe Ree ROD OL PH 


DETERMINES STREAMFUNCTIONS AND GRID POINT VELOCITIES FROM 
POTENTIAL FLOW. 


100 
200 


280 


400 
500 


600 


700 
800 


888 
890 
900 


DIMENSION UA(IXX,1ZZ),VA(IXX,1ZZ),PS(81,52) 
COMMON XP( 10000), ZP( 10000) , JFLAG( 10000), 
1DELX( 10000) ,DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(SO00) ,CT(5000), 


4CMX (5000) ,CMZ(5000) ,XS(10),INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 


NEUSIE TES AW GINES 

DATA PI/3.1415926/ 
IGRIDX=L/INCX+14 
JGRIDX=L/(2*INCX )+1 
IGRIDZ=HPRIME/INCZ+1 
LD2=L/2 

S=2*HPRIME/L 


K=O, DIFFUSION ONLY 


IF(K.EQ.0) GO TO 600 
DOw200) J=1_ IGRIDZ 
ZING Zale 
A=SIN(2*PI*Z/HPRIME ) 
B=SIN(PI*Z/HPRIME ) 

DO 100 1=4, JGRIDX 

SUS TENOR (C1 9] 
PS(I,J)=C*(A*SIN(PI*X/LD2 )+B*SIN(2*P1I*X/LD2) ) 
K=IGRIDX+1-I 
PS(K,JU)=-PS(I,uU) 
CONTINUE 

CONTINUE 
IGRIDX=IGRIDX-1 
JUGRIDZ=IGRIDZ+1 

DO 500 J=1,UGRIDZ 
Z=INCZ*(UJU-1) 
XCPT=LD2-Z/S 

DO 400 1=2,IGRIDX 
X=INCX*( 1-1) 
A=PS(1I,uU-1) 

IF(JU.EQ.1) A=0O 
B=PS(I,uU+1) 

EEGUNEO] UGRIDZ) B=PS(1,J)+(PS(1,0)-PSCL.U=1)) 
UA(I,JU)=(A-B)/(2*INCZ) 
VA(1,J)=(PS(I+1,U)-PS(I-1,U))/(2*INCX) 
CONTINUE 

CONTINUE 

GO TO 800 

DOM/OONL 1A LGRIDx 

DO 700 J=1,IGRIDZ 
UA(I,JU)=0O 

VA(1I,J)=0O 

CONTINUE 

IF(KKK.NE.O) GO TO 900 
IGRIDX=IGRIDX+1 
IGRIDZ=IGRIDZ+1 

DO 890 J=1,IGRIDZ 
WRITE(4,888)(PS(I,U),1=1, 1GRIDX) 
FORMAT(100(1X,F8.5)) 
CONTINUE 

RETURN 

END 
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SUBROUTINE SOURCE(SDEVX,SDEVZ,ITM,IXX,I2ZZ) 


WRITTEN BY: R.RUDOLPH 


DETERMINES AN INITIAL PARTICLE CONFIGURATION, ASSUMED TO BE APPRO- 
XIMATELY GAUSSIAN, FOR EACH POINT SOURCE AT EACH TIMESTEP. 


REAL*8 ITM 

COMMON XP(10000), ZP( 10000) , UFLAG( 10000), 

1DELX( 10000) ,DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(5000) ,CT(5000), 
1CMX (5000) ,CMZ(5000) ,XS(10),INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 
4DELT,T,NS 

DATA IFLAG/O/,PWT/O0.2/ 


CONSTANTS: 

NP= NUMBER OF PARTICLES PRODUCED SO FAR 

XO= X-AXIS SOORCE COORDINATE 

USTAR,STAB ARE FRICTION VELOCITY AND STABILITY (Z/L) 
VK= VON KARMEN COEFFICIENT 

DZ ZS aVERTUCALIEDD YD TEEUST VT TY 

DKXX=HORIZONTAL EDDY DIFFUSIVITY 


ARRAYS: 

Kee Ce eAREe ARTE CEES POS NiPhONS 

VAR IS VARIANCE OF EACH PUFF 

Nite tsSsNUMBERSORS PARTICLES IN EACH PUFF 

NPIT IS NUMBER OF PARTICLES PRODUCED PRIOR TO EACH PUFF 
Ci. lSmiLesCREATION SreiMEs OF. EACH PUBE 


GGNQF IS A SYSTEM SUBROUTINE THAT COMPUTES NORMALLY DISTRIBUTED RANDOM 
NUMBERS, WITH MEAN ZERO AND STANDARD DEVIATION ONE. I.E. N(O,1) 


LD2=L/2 

S=2*HPRIME/L 

DOF 300 K=1,NS 
NT=NT+4 

XO=XS(K) 
ZO=(XO-LD2)*S 
GX Om aie D2) ZO — S20 
SZ=O 

SX=0 

SVAR=O 


DETERMINE NUMBER OF PARTICLES PRODUCED AND THEIR POSITIONS, ASSUMING 
THAT THEY ARE PART OF A NORMAL DISTRIBUTION. 


N=5.24/PWT*DELT*(0.04-0.0075*T/3600) 
ECM LY. 2) Wee 
LL=NP+1 
100 Z=GGNQF (ITM) 
Z=SDEVZ*Z+ZO 
A=LD2-2/S 
X=GGNOF (ITM) 
X=SDEVX*X+XO 
TFC.NOT.(X.GE.A-AND.X.LE.L-A)) GO TO 100 
ELL ea 
ZP(LL)=Z 
SX=SX+X 
SZ=SZ+Z 
IF(LL-NP.GE.N) GO TO 200 
C=C 
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GO TO 100 
Cc 
C CHECK: IS VARIANCE OF SOURCE DISTN. LARGE ENOUGH FOR THE PARTICLES 
C TO BE INCLUDED IN DIFFUSION ROUTINE? IF YES,VAR=-1 
c 
200 CMX(NT)=SX/N 
CMZ(NT)=SZ/N 
NNP=NP+N 
NPP 1=NP+1 
DO 250 I=NPP1,NNP 
SVAR=SVAR+(ZP(I)-CMZ(NT))*(ZP(1I)-CMZ(NT) )+(XP(I)-CMX(NT))* 
1(XP(1)-CMX(NT)) 
250 CONTINUE 
VAR(NT)=SVAR/N 
IF(VAR(NT).NE.O) GO TO 260 
IFLAG=1 
WRITE(7,111) NT,N 
14141 FORMAT(1X,’VARIANCE OF PUFF’,1X,14,1X,’IS ZERO. NUMBER OF PART’, 
HAC RESmT NePUF Re Sed Xen ri o)) 
260 SIG=SORT(VAR(NT)) 
AA=AMAXO(INCX,INCZ) 
IF(SIG.GT.AA) VAR(NT)=-1 
NIT(NT)=N 
NPIT(NT)=NP 
NP=NP+N 
CmaNin aii 
WRITE(6,333) NT,T,VAR(NT),NPIT(NT),NIT(NT),NP 
333 FORMAT(3X,14,2(3X,F8.1),3(3X,15)) 
IF(IFLAG.NE.1) GO TO 300 
WRITE(8,222)T,NT,K 
222 FORMAT(//3X,’PREMATURE TERMINATION IN SUBROUTINE SOURCE. CHECK 
(ERROR FLU ERadXe nT IME=" 4X oh@n 11x, 2 FOR PURE 1xe 158 4X%e" FORESOU 4 
OERGEMEOGAI LON Gat xl 2/7) 
STOP 
300 CONTINUE 
RETURN 
END 
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SUBROUTINE DIF 1(1II,DKMX,DKMZ,IXX,1ZZ,UD,VD,CONC) 


WRITTEN BY: R.RUDOLPH 


COMPUTES DIFFUSION VELOCITIES AT GRIDPOINTS 


0000 0 


DIMENSION UD(IXX,1Z2Z),VD(1IXX,1ZZ),CONC(IXX,1ZZ) 

COMMON XP( 10000) ,ZP( 10000) , JFLAG( 10000) , 
2DELX( 10000) ,DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(5000) ,CT(5000), 
3CMX (5000) ,CMZ(5000) ,XS(10),INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 
4DELT,T,NS 

DATA VK/.35/,AA/1.0/ 

LOA=L /2 

IGRIDX=L/INCX 

IGRIDZ=HPRIME/INCZ+1 

S=HPRIME*2/L 


@ 
C 1S DIFFUSION TO BE NEGLECTED (EXCEPT INITIALLY) ? 
(e; 


IF(AA.GE.1.0) GO TO 30 
DO 20 J=1,IGRIDZ 
DONZO. 1=41 TGRIDX 
UD(1I,JU)=0 
WB 1 oS NEO 
20 CONTINUE 
GO TO 800 
30 B=1+4.7*STAB 
TROSHABLLT-O) B=(1—15*STAB)+*(—0025) 
IFLAG=O 
DO 500 J=1,IGRIDZ 
Z=INCZ*(uU-1) 
DO 460 I=2,I1GRIDX 
X=INCX*( 1-1) 
IF(III.EQ.0) GO TO 40 
SZ=(LD2-X)*S 
ie (OS, CF LD) Svze-Sz 
SZ=ABS(Z-SZ) 
DKZZ=VK*USTAR*SZ/B 
DKXX=DKZZ*SZ**(0.07) 
GO TO 50 
40 DKXX=DKMX 
DKZZ=DKMZ 
50 VM=ABS(2*DKZZ/INCZ) 
UM=ABS(2*DKXX/INCX ) 
A1=CONC(I,JU) 
A2=CONC(1-1,U) 
IF(JU.NE.IGRIDZ) GO TO 60 
MIECONEC I, GW) 
A2=CONC(I-1,U-1) 
60 A@=CONC(1I-1,U-1) 
A4=CONC(1,U-1) 
He(d NE.) Co WO 7) 
A3=CONC(I-1,U) 
A4=CONC(I,U) 
70 IF((Z-INCZ).GE.HPRIME) GO TO 100 
Cc 
C FOR INTERSECTION BELOW VALLEY BOUNDARY 
G 
Z1=(X-LD2)*S 
1EOC Le De) Ze S744 
ie (2 =—72)) Le ©) Cle We) Wee 
UD(1I,J)=0O 
VD(I,JU)=0O 
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GO TO 460 
Cc 
C CHECK - ALL CONCENTRATIONS SHOULD BE POSITIVE 
Cc 


VO) We Wo. CAD ir ©. 0. Ae. LPO. OR Ne). Li ©. O22 har ©) Co WO Bee 
Wile E(i7entid in) elec vA AGM AQ mIAIA 
111  FORMAT(3X,’ERROR IN DETERMINING GRIDPOINT DIFFUSION VELOCITY 
4 FOR GRID’ ,2X,13,’,’,13,2X,’- NEGATIVE CONCENTRATION ENCOUNTERED’ // 
ON an Smon Nl er Qn oa AS =n CBE Gro Ad = 78 F Sia) 7s) 
UD(1,U)=-999999 
VD(I,JU)=-999999 


IFLAG=1 

GO TO 460 
Cc 
CEGARCUEAT ESD LEEUSTONS VELOCT TlrESm(REREE LANGER Oise eS)} 
Cc 


200 D=A1+A2+A3+A4 
IF(D.EQ.0) GO TO 400 
UD(I,JU)=2*DKXX* (A2+A3-A1-A4)/(INCX*D) 
VD(I,U)=2*DKZZ*(A3+A4-A1-A2)/(INCZ*D) 
GO TO 460 
400 UD(I,uJ)=0 
VOL Ee 
460 CONTINUE 
500 CONTINUE 
IF(IFLAG.NE.1) GO TO 800 
WRITE(8.777) T,NT 
777. FORMAT(//1X, ‘PREMATURE TERMINATION IN DIF1‘/10X/‘AT TIME’, 
fine Foe. “FOR PUFF’ 41X,17. CHECK ERROR FILE’//) 
STOP 
800 RETURN 
END 
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SUBROUTINE SUM(UJUJU,%,Z,1XX,1ZZ,CONC, AREA) 
C WRITTEN BY:R.RUDOLPH 
Cc 


C DETERMINES CELL PARTICLE CONCENTRATION 
Cc 


DIMENSION CONC(IXX,1ZZ),AREA(IXX,1ZZ),1ISUM(81,52) 


COMMON XP( 10000), ZP( 10000) , UJFLAG( 10000), 


1DELX (10000) ,DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(5000) ,CT(5000), 
4CMX (5000) ,CMZ(5000) ,XS(10),INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 


1DELT,T,NS 
II=X/INCX+1 
JJ=Z/INCZ+1 
TJU=(X+150)/INCX+14 
JI=(X-250)/(8*INCZ)+14 


TF ONT.EQs 1) MP WRITE(G)300) AREA(II, JU) ,AREACTU, JI) 


IGRIDX=L/INCX 
IGRIDZ=H/INCZ 
C 
C INITIALIZE PARTICLE/CELL ARRAY 
c 
DONQ4ONl= IUGR IDX 
DO 240 J=1,IGRIDZ 
ISUM(I,JU)=0O 
240 CONTINUE 
DO 250 I=1,NP 
INT=XP(1)/INCX+1 
JUNT=ZP(1)/INCZ+14 
ISUM(INT,UNT)=ISUMC INT, UJNT)+1 
250 CONTINUE 
G 
C CALCULATE CELL PARTICLE CONCENTRATION 
G 
DO 290 I=1,IGRIDX 
DO 290 J=1,IGRIDZ 
IE(AREA(I,J).LT.1.0E-4) GO TO 260 
CONC(I,JU)=ISUM(1I,JU)/AREA(I,U) 
GO TO 290 
260 CONC(I,J)=-1 
290 CONTINUE 
le (dWh eo) Ee ie) 2dse@ 
IITI=MOD(NT, JuJu) 
me(iits .Ne.@)) Co wa see 


WRITE(3,300)CONC(II, UJ) ,CONC(II,JUJU+1),CONC(II,JJ+2), 


4CONC(IU,JI),CONC(IU,JUI+1),CONC(IJU,JI+2),T 
300 FORMAT(1X,6F9.5,1X,F8.1) 
350 RETURN 

END 
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SUBROUTINE INTER(SDEVX,SDEVZ,III,DKMX,DKMZ,INV,IXX,1ZZ,UD, 
1VD,UA,VA) 


WRITTEN BY: R.RUDOLPH 


INTERPOLATES VELOCITIES FROM GRID PTS.TO PARTICLE POSITIONS 
DIFFERENTIATING BETWEEN "NORMAL" AND SUBGRID SCALE DIFFUSION 
CALCULATES MAX. GRID VELOCITY AND LENGTH OF SUBSEQUENT TIMESTEP 


SYMBOLS: 

UA,UD : ADVECTION AND DIFFUSION VELOCITIES INTERPOLATED FROM GRID 
INTERSECTIONS 

PHIM : DIMENSIONLESS WIND SHEAR 

ERO@eRAl EOD 1 SSiPAT TON: OF SEDDY ENERGY 

VK : VON KARMAN COEFFICIENT 


DIMENSION UD(IXX,1ZZ),VD(1XX,1ZZ),UA(IXX,1ZZ),VA(IXX,1ZZ) 
COMMON XP( 10000), ZP( 10000) , JFLAG( 10000), 

1{DELX( 10000) , DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(5000) ,CT(5000), 
1CMX (5000) ,CMZ(5000) ,XS(10),INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 
4DELT,T,NS 

REAL UDX,UDZ,UAX,UAZ 

DATA VK/.35/,BB/1.0E-4/,CC/1.0E-7/,EE/O. 10/ 

UMAX=0O 

VMAX =O 

AMUAZ=O 

IFLAG=O 

LOeG=L/2 

S=HPRIME/LD2 

PHIM=1+4.7*STAB 

IF(STAB.LT.O) PHIM=(1-15*STAB)**(-.25) 

DO) 500 1=1,NT 

LLeNey EW, 

DOM OOM =a lal: 

II=I+NPIT(L1) 

INT=XP(II)/INCX+1 

JNT=ZP(11)/INCZ+1 

X1=XP(I1)-(CINT-1)*INCX 

X2=INCX-X1 

Z1=ZP(1I1)-(JUNT-1)*INCZ 

ZI NCZAZ A 


INTERPOLATE ADVECTION VELOCITIES TO PARTICLE POSITION 


UAX=Z2* (UAC INT, UNT ) *X2+UA( INT+1,UNT )*X1) 

UAX =UAX+Z1* (UAC INT, UNT+1)*X2+UA(INT+1,UNT+4)*X1) 
UAX=UAX/( INCX*INCZ) 

UAZ=X2* (VAC INT, UNT )*Z2+VAC INT, UNT+1)*Z1) 
UAZ=UAZ+X1* (VAC INT+1, UNT+1)*Z1+VACINT+14,UNT)*Z2) 
UAZ=UAZ/(INCX*INCZ) 

A=ABS(UAZ) 

AMUAZ=AMAX 1(AMUAZ,A) 


TEST VARIANCE : IF POSITIVE USE SUB-GRID DIFFUSION 
IF(VAR(L1).GE.0) GO TO 200 

INTERPOLATE DIFFUSION VELOCITIES TO PARTICLE POSITION 
UDX=Z2* (UD( INT, UNT ) *X2+UD(INT+1,UNT )*X1) 
UDX=UDX+Z1* (UD( INT, UNT+ 1) *X2+UD(INT+1,UNT+1)*X1) 


UDX=UDX/ ( INCX*INCZ) 
UDZ=X2*(VD( INT, UNT )*Z2+VD( INT, UNT+1)*Z1) 
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UDZ=UDZ+X1* (VDC INT+1,JUNT+1)*Z1+VD(INT+1,UNT)*Z2) 
UDZ=UDZ/( INCX*INCZ) 

DELX(1I1I)=(UDX+UAX ) *DELT 

DELZ(I1I)=(UDZ+UAZ)*DELT 

GORTORSOO 


ASSUME PUFF TO BE IN CONSTANT (SURFACE) STRESS LAYER,AND TURBULENCE 
TO BE IN THE INERTIAL SUBRANGE 


CALCULATE THE DISTANCE DIFFUSED BY THE PARTICLE 
RE ANGE a Simro © 


OOOO OOS 


XOXO) WE WI. (2@),C)) EO wel Behl 
Z=(LD2-XP(I1))*S 
miei), Gir Loe@)) 74ss72 
Z=Z2P\C1 W) =z 
TR@ZASGE +O) GORTON240 
WRIMECGim222)) =xPiChID ZeChi eZ UA UAZ 
BAD JOM VNU Ck 24s AR? 9 I, Bie teh, We, 274 — CR TO eh, ha UU 6g Pe Ie ey) 
STOP 
240 IF(ABS(Z).LT.BB) GO TO 250 
EPS=USTAR*USTAR*USTAR/(VK*Z)*(PHIM-STAB ) 
D=T-CT(L1)+1.5*(SDEVZ*SDEVZ/EPS) **(.333) 
DD=T-CT(L1)+1.5*(SDEVX*SDEVX/EPS)**(.333) 
ee) 10 2O 
D5 OmD =~ Cir Clalb) 
DD=D 
276) WE NOT. (Li .Ce.02. Lr Kee) ele) We So) 
DELZ( II )=UAZ*DELT 
DELX(II)=UAX*DELT 
GO TO 300 
290 DELZ(II)=UAZ*DELT+(((1+DELT/D)**(1.5))-1)*(ZP(1I1I)-CMZ(L1)) 
DELX(II)=UAX*DELT+(((1+DELT/DD)**(1.5))-1)*(XP(1II)-CMX(L1)) 
GO TO 300 
291 DX=SDEVX*SDEVX/(2*DKMX )+T-CT(L1) 
IF(DX.EQ.0) GO TO 292 
DELX( II )=UAX*DELT+(XP(II)-CMX(L1))*(SQRT(1+DELT/DX )-1) 
GO 10) 294 
DEQ INZY 
B=XP(1I1)-CMX(L1) 
IF(B.EQ@.0) GO TO 293 
DELX(II )=INCX/2*SIGN(A,B) 
GOmnom294 
293 DELX(II)=UAX*DELT 
294 DZ=SDEVZ*SDEVZ/(2*DKMZ)+T-CT(L1) 
IF(DZ.EQ.0) GO TO 295 
DELZ(II )=UAZ*DELT+(ZP(11)-CMZ(L1))*(SQRT(1+DELT/DZ)-1) 
GO TO 300 
295 A=1 
B=ZP(1I1)-CMZ(L1) 
IF(B.EQ.0) GO TO 296 
DELZ(II)=INCZ/2*SIGN(A,B) 
Go 18 see 
296 DELZ(II)=UAZ*DELT 
c 
C FIND MAXIMUM VELOCITY ON THE GRID 
G 
300 U=ABS(DELX(II)/DELT) 
V=ABS(DELZ(I1)/DELT) 
UMAX=AMAX 1(UMAX,U) 
VMAX=AMAX 1(VMAX,V) 
400 CONTINUE 
500 CONTINUE 


Cc 
C DETERMINE THE LENGTH OF THE NEXT TIMESTEP 
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C AND THE RATIO OF DIFFUSION TO TOTAL VELOCITIES 


(e 
IF(VMAX.EQ.0) VMAX=1.0E-5 
V=(VMAX-AMUAZ ) / VMAX 
TE (VEO. O) V=1.0E-5 
INV=1.0/V 

Cc 


A1=INCX/(2*UMAX ) 
A2=INCZ/(2*VMAX ) 
DT=AMINi(A1,A2) 
A3=DELT*1.3 
DELT=AMIN1(DT,A3) 
IF(DELT.LT.EE) IFLAG=1 
WRITE(6,666)A3,DT,UMAX, VMAX ,DELT 
666 FORMAT(5(3X,F10.3)) 
ie (WPLAKE ME.) GB) he) Wee) 
WRITE(8,777) T,EE 
777 FORMAT(//1X, ‘PROGRAM ABORT IN SR INTER AT TIME’,1X,F8.1,5X, 
4’/TIMESTEP LESS THAN ‘,1X,F8.5//) 
STOP 
700 RETURN 
END 
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SUBROUTINE TSTEP(IOUT,NOTIN) 


C WRITTEN BY:R.RUDOLPH 


G 


C ADVANCES PARTICLES ALONG STREAMLINES, REFLECTS PARTICLES 


Cc 
Cc 


C 


AT 


APPROPRIATE BOUNDARIES, AND RECALCULATES THE VARIANCE OF EACH PUFF. 


COMMON XP( 10000) ,ZP( 10000) , UFLAG( 10000), 

IDELX( 10000) , DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(5000) ,CT(5000), 
1CMX (SOOO) , CMZ(S000) ,XS( 10), INCX, INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 
1DELT,T,NS 

DATA IFLAG/O/,PI/3.1415926/,BB/1.00E-4/,CC/1.00E-3/ 

LDQ hyo 

S=2*HPRIME/L 

BO=ARSIN(S) 

DOM COOn =m Nil 

eo 

SZ=© 

I=IN (Le). (20) @))) Cle) lel See 

LIL SNL 3h) 

DORSCOn lA Le 

II=I+NPIT(L1) 

HELA it NEO 

XP1=XP(II1) 

ZP1=ZP(II1) 


C DETERMINE NEW PARTICLE POSITIONS 


Cc 


CGI (@) 


C 


XP2=XP(11)+DELX(II) 
Ape 7 eiate) DE eZ (iet)) 


ChECkeas IS) PARTICEE OUT TOR BOUNDS? ~IF*VES*PUFEAG=4 


10 


IF((XP2.LT.INCX).OR.(XP2.GT.(L-INCX))) JFLAG(II)=14 
WE (7222) CE. ©.) EO We Ke 

JFLAG(II)=1 

GO 70: 30 

IF(ZP2.LE.HPRIME) GO TO 90 


C PARTICLE TRANSPORTED ABOVE RIDGELINE 


C 


C 


20 


30 


sis) 


NOTIN=NOTIN+1 
we (NOTIN. Lip. KOUn)) Ca) yo 2o 

UF EAG(IT)=4 

NOTIN=O 

GO TO 80 

ZP2=2*HPRIME-ZP2 

IF(ZP2.GT.HPRIME) UFLAG(II)=1 
IF(ZP2.EQ.ZP1) JFLAG(II)=1 
IF(JUFLAG(II).EQ.1) GO TO 95 
X1=LD2-ZP2/5S 

THGCeom hex deOREXP2) Gil —41)) GO miOn 100 
XP2H=XP2 

ZP2H=ZP2 

CORIO MZOC 


C PARTICLE ADVECTED INTO GROUND. FIND GROUND-TRAVECTORY INTERCEPT 
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100 


LECKPA SEG XP2)EGD 7O AIO 
S$1=(ZP1-ZP2)/(XP1-XxP2) 
S2=S 

LE(Xe27U1.LD2) S2=-s 
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IF S1=S2, PARTICLE IS AT GROUNDLINE, AND THEREFORE IN THE ATMOSPHERE 


ECS. 2@).S2)) Ee) ie) Fore) 
B2=HPRIME 
IF(S2.GT.0) B2=-B2 
X=(S1*XP14+B2-ZP1)/(S1-S2) 
Z=(S2*ZP1-S1*B2-S1*S2*XP1)/(S2-S1) 
GO TO 140 
1Omxe=xXP2. 
Z=(X-LD2)*S 
PRO LP (BE) yea sre 
4 © mein (Gen Gila ae AN Din xen Gil Xe 2.) Oke Oxon ate PS ANDE Xe lex e2) ORE 
GAR Giz PAP ANDEZ  GiazZe2)mORMCZ. Ll eZee ANDe Zt Ze2)) GOn 10M 145 
GO TO 160 
IAS MRIECY 5 VSO MEW 5 7424) Be. AA) 9 72 
1450 FORMAT(3X, ‘UNREASONABLE INTERSECTION’ ,3X, ‘BEGIN=’,E23.16, 7’, 
Eo Smo SxXaneEND=4nE 20m 1Gmeu FOS 81Gs3aXs4INTERSEGH=", E2016. 2 
2E23.16) 
IF (ABS(XP1-X).GE.CC.AND.ABS(XP2-X).GE.CC.AND.ABS(ZP1-Z).GE.CC 
1.AND.ABS(ZP2-Z).GE.CC) GO TO 155 
JFLAG(II)=1 
XP2H=XP2 
ZP2H=ZP2 
GO TO 700 
155 IFLAG=1 
(fo) (0) TKOO 


REE ERECT AIRES PAR i Cle 


160 A1=COS(BO) 
NDES 
D1=SORT((XP1-XP2)*(XP1-XP2)+(ZP1-ZP2)*(ZP1-ZP2) ) 


IF D1=O, PARTICLE IS ON THE GROUND 


i= ((13)9 .{2@)@)) Ea} syle) 7hoke 
D3=SORT( (X-XP2)*(X-XP2)+(Z-ZP2)*(Z-ZP2)) 
TECXP 2G O2) GOMICE200 
X=LD2-X 
XP 1=LD2-XP1 
XP2=LD2-XP2 
Why Sz 
GO TO 300 
DOO WE XK= Xe ED2 
XP1=XP1-LD2 
XP2=XP2-LD2 
JJ=O0 


ROTATE AXIS BY BO RADIANS 


300 XP1R=XP1*A1+ZP1*A2 
XP2R=XP2*A1+ZP2*A2 
X1=X*A1+Z*A2 
Z1=-X*A2+Z*A1 

MEGZA ee OMe Zall= OF 
ARG=(XP1R-XP2R)/D1 
IF(ABS(ARG).GT.1.0) ARG=SIGN(1.0,ARG) 
B=ARCOS(ARG) 

IF(XP1R.LT.XP2R) GO TO 400 

I1=-1 

GO TO 500 

400 I11=1 
B=PI-B 

500 XP2R=X1+11*D3*COS(B) 
ZP2R=Z1+D3*SIN(B) 
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C ROTATE AXIS TO ORIGINAL POSITION 
G 
XPN=XP2R*A1-ZP2R*A2 
ZPN=XP2R*A2+ZP2R*A 1 
UF dh) GE. 1) GO Fa Geo 
XPN=XPN+LD2 
GO TO 650 
600 XPN=LD2-XPN 
650 XP2H=XPN 
ZP2H=ZPN 
A2=LD2-ZP2H/S 
IF (XP2H.GE.A2.AND.XP2H.LE.L-A2) GO TO 700 
TECXP2HeGieLD2)) GOTO 66o 
R=(XP2H-A2)/(XP2H+1) 
GO TO 670 
660 R=(XP2H-L+A2)/XP2H 
670 IF(ABS(R).GT.BB) GO TO 680 
UREAGGDI=x 
GO TO 700 
680 IFLAG=1 
Wee fixe 1a ZPiaXk, 2a XP 2. 2P2 XPAR XP2R % 1.7215 XP2H ZPOH 
Ti TREORMAT MIXES REF L BELOW GND’, 1X.?XP4,ZP4". 1%, 2-12.38) 1x, 
tree 2 ioe Oya XPD. ePO”. 4X OF 19 Be AX Pe XPA pop 4p) 
De eo eX eA AX OF ADB / 1k,  XPIH) TPO Ce 
32F12.8) 
700 XP (Il )=XP2H 
ZP(II)=ZP2H 
SX=SX+XP2 
SZ=SZ+ZP2 
800 CONTINUE 
Cc 
C IF DISTRIBUTION IS SUB-GRID SCALE, CALCULATE CMX,CMZ 
C CMX,CMZ ARE COORDINATES OF THE CENTER OF MASS OF THE PUFF 
Cc 
IF(VAR(L1).LT.O) GO TO 950 
GUL VES EL 
CMZ(L1)=SZ/LL 
C CALCULATE NEW PUFF VARIANCE 
SVAR=O 
850 DO 900 I=1,LL 
II=I+NPIT(L1) 
SVAR=SVAR+(XP( II )-CMX(L1))*(XP(II)-CMX(L1))+(ZP(II)-CMZ(L14)) 
1 CZPiGL TH) —CMZ (1 ))) 
800 CONTINUE 
SVAR=SVAR/LL 
AA=AMAXO(INCX, INCZ) 
IF(SQRT(SVAR).GT.AA.OR.SQRT(SVAR).LT.BB) SVAR=-4 
VAR(L1)=SVAR 
GO TO 9g60 
950 VAR(L1)=-1 
960 IF(IFLAG.NE.1) GO TO 1000 
WRITE(8,999) NT,T 


999 FORMAT(//3X, ‘PREMATURE TERMINATION IN SUBROUTINE TSTEP AT TIMESTEP 


I? 9 Wk, MSO, LIND ite”, i RS}. 1) SM COREG faleletole (rullie? 7/7) 
STOP 
1000 CONTINUE 
RETURN 
END 
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SUBROUTINE RENUM 


WRITTEN BY: R.RUDOLPH 


RENUMBERS EXISTING PARTICLES, DESTROYING THOSE FOUND TO BE OUTSIDE 
THE GRID SYSTEM ( VIA JFLAG ARRAY ) 


COMMON XP( 10000) ,ZP( 10000) , JFLAG( 10000), 

1DELX (10000) ,DELZ( 10000) , VAR(5000) ,NIT(5000) ,NPIT(5000),CT(5000), 

1CMX (5000) ,CMZ(5000) ,XS(10),INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 

1DELT,T,NS 

DIMENSION TX( 10000) ,TZ( 10000) ,DF 1(5000) 

N=0O 

JSq 

DF1i(1)=0 

DOMZOOn EH = IeNit 

IF(NIT(L1).EQ.0) GO TO 150 

fale Neate) 

OY) Ge) re) (LL 

II=I+NPIT(L1) 

IF(JFLAG(II).EQ.1) GO TO 50 

PMC UNE MeC TT”) 

TZ(Y Y=EZP(IL) 

J=J+14 

GOMTGMOO 

NIT(L1)=NIT(L1)-1 

N=N+4 

WRITE(6,111) XP(II),ZP(II) 

FORMAT(41X, ‘DETECTED OUT-OF-BOUNDS X=’,1X,F8.4,2X,’Z=’,1X,F8.4) 
CONTINUE 

DF1(L1+1)=NPIT(L1+1)-N 

CO) TO Boo 

DF1(L1+1)=DF1(L1) 

CONTINUE 

NP=uJ-1 

DOMSOOm es = iieNil 

NPIT(L1)=DF1(L1) 
WRIGEC CRASS IRIE CEC VAR COI) RCMx CLI) RCMZ (ld a NPIL CEI c 
4NIT(L1),NP 
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“500 


600 
700 


CONTINUE 

IF(NP.LT.1) GO TO 700 
DO 600 I=1,NP 

OCT NST IE) 
ZP(1)=TZ(1) 

CONT INUE 

RETURN 

END 
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SUBROUTINE OUTPT(IXX,1ZZ,CONC) 


WRITTEN BY: R. RUDOLPH 


WRITES CONCENTRATION AND PARTICLE POSITION DATA INTO FILES (12,13) 
FOR USE IN PLOTTING PROGRAMS. 


DIMENSION CONC(IXX,IZZ) 
COMMON XP( 10000) ,ZP( 10000) , JFLAG( 10000), 
1DELX( 10000) ,DELZ( 10000) , VAR(5000) ,NIT(S5000) ,NPIT(5000) ,CT(5000), 


1CMX (5000) ,CMZ(5000) ,XS(10),INCX,INCZ,HPRIME,H,L,NP,USTAR,STAB,NT, 


1DELT,T,NS 

DATA ZAP/‘***// 
GRD X= l/ALNGX 
IGRIDZ=H/INCZ 


WRITE(12,900)NP 
FORMAT(I5) 


WRITE(12)(XP(1I),ZP(1),1=1,NP) 
WE HE Va, Waal) 727M INT 5 TF 


DO 1400 J=1,IGRIDZ 

WRITE(13, 1300) (CONC(I,JU),1=1,IGRIDX) 
FORMAT(80(1X,F8.4)) 

CONTINUE 


WRITE( 13,1444) ZAP,NT,T 

FORMAT(1X.A4, 1X,’ TIMESTEPS=", 1X,15,3X, ’TIME=”, 1X,F8.1) 
RETURN 

END 
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